{"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_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Imports","metadata":{}},{"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","metadata":{"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import os\nimport cv2\nimport pydicom\nimport pandas as pd\nimport numpy as np \nimport tensorflow as tf \nimport matplotlib.pyplot as plt \nimport random\nfrom tqdm.notebook import tqdm \nfrom sklearn.model_selection import train_test_split, KFold\nfrom sklearn.metrics import mean_absolute_error\nfrom tensorflow_addons.optimizers import RectifiedAdam\nfrom tensorflow.keras import Model\nimport tensorflow.keras.backend as K\nimport tensorflow.keras.layers as L\nimport tensorflow.keras.models as M\nfrom tensorflow.keras.optimizers import Nadam\nimport seaborn as sns\nfrom PIL import Image\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    \nseed_everything(42)","metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"config = tf.compat.v1.ConfigProto()\nconfig.gpu_options.allow_growth = True\nsession = tf.compat.v1.Session(config=config)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/train.csv') ","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 1. Background: general information about your chosen ADS\n\n* **What is the purpose of this ADS? What are its stated goals?**\n\nThe purpose of this ADS is to automatically predict a patient’s severity of decline in lung function based on a CT scan of their lungs. ","metadata":{}},{"cell_type":"markdown","source":"# 2. Input and output\n* **Describe the data used by this ADS. How was this data collected or selected?**\n\nFor each patient, there is a baseline chest CT scan and clinical information. The baseline chest CT scan is at Week 0 and over 1-2 years, there are followup visits with spirometry tests. There are train and test CSV files, which hold the unique patient ID, the number of weeks before or after the CT scan for a visit where FVC value (mL) was determined, “a computed field which approximates the patient’s FVC as a percent of the typical FVC for a person of similar characteristics,” age, sex and smoking status. The details to the collection is not explicitly stated. The only detail mentioned about the data is that it is real medical data.","metadata":{}},{"cell_type":"markdown","source":"* **For each input feature**\n    * Describe its datatype\n    * Give information on missing values\n    * Value distribution\n    * Show pairwise correlations between features if appropriate. \n    * Run any other reasonable profiling of the input that you find interesting and appropriate.\n    * What is the output of the system (e.g., is it a class label, a score, a probability, or some other type of output), and how do we interpret it?","metadata":{}},{"cell_type":"code","source":"# Data types\n\ntrain.info()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Missing Values\n\ntrain.isna().sum()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Value distribution\n\nplt.figure(figsize=(20,15))\n\nplt.subplot(2, 3, 1)\nplt.hist(train['Weeks'])\nplt.title('Week')\nplt.xlabel('Week Value')\nplt.ylabel('Count')\n\nplt.subplot(2, 3, 2)\nplt.hist(train['FVC'])\nplt.title('FVC')\nplt.xlabel('FVC Value')\nplt.ylabel('Count')\n\nplt.subplot(2, 3, 3)\nplt.hist(train['Percent'])\nplt.title('Percent')\nplt.xlabel('Percent Value')\nplt.ylabel('Count')\n\nplt.subplot(2, 3, 4)\nplt.hist(train['Age'])\nplt.title('Age')\nplt.xlabel('Age Value')\nplt.ylabel('Count')\n\nplt.subplot(2, 3, 5)\nplt.hist(train['Sex'])\nplt.title('Sex')\nplt.xlabel('Sex Value')\nplt.ylabel('Count')\n\nplt.subplot(2, 3, 6)\nplt.hist(train['SmokingStatus'])\nplt.title('Smoking Status')\nplt.xlabel('Smoking Status Value')\nplt.ylabel('Count')\n\n\nplt.show()\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Pairwise Correlations\n\nplt.figure(figsize=(12,10))\ncor = train.corr()\nsns.heatmap(cor, annot=True, cmap=plt.cm.Reds)\nplt.title('Correlation Matrix - Train DF')\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"imdir = \"/kaggle/input/osic-pulmonary-fibrosis-progression/train/ID00123637202217151272140\"\nprint(\"total images for patient ID00123637202217151272140: \", len(os.listdir(imdir)))\n\n# view first (columns*rows) images in order\nfig=plt.figure(figsize=(12, 12))\ncolumns = 4\nrows = 5\nimglist = os.listdir(imdir)\nfor i in range(1, columns*rows +1):\n    filename = imdir + \"/\" + str(i) + \".dcm\"\n    ds = pydicom.dcmread(filename)\n    fig.add_subplot(rows, columns, i)\n    plt.imshow(ds.pixel_array, cmap='gray')\nplt.title('Example CT Scan Images for a Patient')\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 3. Implementation and validation:\n* **Describe data cleaning and any other pre-processing**\n\nThe only data cleaning done for this project is omitting two patient ids and their corresponding images. Some of the image files associated with Patient IDs 'ID00011637202177653955184' and 'ID00052637202186188008618' are corrupted and cannot be loaded. \n\nThe pre-processing is typical for TensorFlow. First, we are given a test df with values corresponding to patients. This is vectorized so it can be used as input to the neural net. The other input to the neural net are the image files, which are 512 x 512 pixels in size. \n\n* **Give high-level information about the implementation of the system**\n\nThe solution in this notebook uses EfficientNet B5 and quantile regression neural networks (QRNN). EfficentNets are a family of CNN models that shares the same convolution operations as the baseline network but the depth, width and resolution are all uniformly scaled with a compound coefficient depending on how much resources are available. It uses significantly less parameters and computing power than comparable models. Quantile regression disregards the assumption of constant vari- ance for the error term typical to linear regression and it will essentially try to fit a line that splits the data such that there is an certain amount of data below that line defined by the quantile (e.g., 0.75); this concept is applied to the loss function for a neural network. EfficientNets utilized CT scans and tabular data, whereas QRNN relied on tabular data. The decisions from these models were blended together as an ensemble to make the final predictions.\n\n* **How was the ADS validated? How do we know that it meets its stated goal(s)?**\n\nThe evaluation metric is based on a modified version of Laplace Log Likelihood and its purpose is to evaluate how far the predicted FVC value is from the true value and the model’s confidence in that prediction. If the predicted value is far from the true value (clipped at 1000) but less confident (it is represented as a high value and can also be interpreted as standard deviation or uncertainty), then the model is not penalized as much. The confidence value is clipped at 70 as that reflects the ”approximate measurement uncertainty in FVC.” If the absolute error is high and the model is very confident, then the model will be penalized harshly.\n","metadata":{}},{"cell_type":"markdown","source":"# Outcomes\n\n* **Analyze the effectiveness (accuracy) of the ADS by comparing its performance across different subpopulations.**\n\n* **Select one or several fairness or diversity measures, justify your choice of these measures for the ADS in question, and quantify the fairness or diversity of this ADS.**\n\n* **Develop additional methods for analyzing ADS performance: think about stability, robustness, performance on difficult or otherwise important examples (in the style of LIME), or any other property that you believe is important to check for this ADS.**","metadata":{}},{"cell_type":"markdown","source":"# Summary\n* **Do you believe that the data was appropriate for this ADS?**\n\nYes, we do believe the type of data is appropriate becasue in the medical field, doctors would be given the same data. CT scans and some historical data. The goal of this ADS is to have a model thats trained on much more data than a doctor can reference by themselves. However, it seems inappropriate that data is 80% male and 20% female. Also, the sample test output only has males, which is also inappropriate.\n\n* **Do you believe the implementation is robust, accurate, and fair? Discuss your choice of accuracy and fairness measures, and explain which stakeholders may find these measures appropriate.**\n\n\n* **Would you be comfortable deploying this ADS in the public sector, or in the industry? Why so or why not?**\n\n\n* **What improvements do you recommend to the data collection, processing, or analysis methodology?**","metadata":{}},{"cell_type":"markdown","source":"# Start of Notebook","metadata":{}},{"cell_type":"code","source":"def get_img(path):\n    d = pydicom.dcmread(path)\n    return cv2.resize(d.pixel_array / 2**11, (512, 512))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Example Image\n# All images are 512 x 512 images\n# The directory structure is that each unique id from train df\n# has a folder and in that folder are different number of image files\n# Some patients had more CT scans and other patients have fewer CT scans\n\nex_img = get_img(f'../input/osic-pulmonary-fibrosis-progression/train/ID00007637202177411956430/4.dcm')\n\nex_img.shape, np.min(ex_img), np.max(ex_img)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Linear Decay (based on EfficientNets)","metadata":{}},{"cell_type":"code","source":"def get_tab(df):\n    vector = [(df.Age.values[0] - 30) / 30] \n    \n    if df.Sex.values[0] == '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    return np.array(vector) ","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"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)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## CNN for coeff prediction","metadata":{}},{"cell_type":"code","source":"def get_img(path):\n    d = pydicom.dcmread(path)\n    return cv2.resize(d.pixel_array / 2**11, (512, 512))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"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(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","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"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):\n    inp = Input(shape=shape)\n    base = get_efficientnet(model_class, shape)\n    x = base(inp)\n    x = GlobalAveragePooling2D()(x)\n    inp2 = Input(shape=(4,))\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#     print(x)\n    model = Model([inp, inp2] , x)\n    \n    weights = [w for w in os.listdir('../input/osic-model-weights') if model_class in w][0]\n#     model.load_weights('../input/osic-model-weights/' + weights)\n    model.load_weights('../input/effnet-b5-30epochs-1/effnet_30.h5')\n    return model\n\nmodel_classes = ['b5'] #['b0','b1','b2','b3',b4','b5','b6','b7']\nmodels = [build_model(shape=(512, 512, 1), model_class=m) for m in model_classes]\nprint('Number of models: ' + str(len(models)))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"models[0].output","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from sklearn.model_selection import train_test_split \n\ntr_p, vl_p = train_test_split(P, \n                              shuffle=True, \n                              train_size= 1) ","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sns.distplot(list(A.values()));","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def 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)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Need to figure out how to do prediction with subset of data","metadata":{}},{"cell_type":"markdown","source":"## Split data based on sex to run algorithm on sub populations","metadata":{}},{"cell_type":"code","source":"# Final DFs for subpopulations\n# Final female and male dfs that are the third param to effnet_iter\n\nBAD_ID = ['ID00011637202177653955184', 'ID00052637202186188008618']\n\nfemale_df = train[train.Sex=='Female']\nfemale_df = female_df[~female_df['Patient'].isin(BAD_ID)]\nfemale_df = female_df[~(female_df.duplicated(['Patient']))].sample(n=37)\n\nmale_df = train[train.Sex=='Male']\nmale_df = male_df[~male_df['Patient'].isin(BAD_ID)]\nmale_df = male_df[~(male_df.duplicated(['Patient']))].sample(n=37)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def make_sub_df(df):\n    tempy = pd.DataFrame(columns=['Patient', 'Weeks'])\n\n    for p in tqdm(df.Patient.unique()):\n        x = [] \n        tab = [] \n        ldir = os.listdir(f'../input/osic-pulmonary-fibrosis-progression/train/{p}/')\n        for i in ldir:\n            tempy = tempy.append({'Patient' : p, 'Weeks': str(i.split('.')[0])}, ignore_index=True)\n            \n    return tempy.sort_values(by=['Patient', 'Weeks'], ascending=[True, True])\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tempy_female = make_sub_df(female_df.copy())\ntempy_male = make_sub_df(male_df.copy())","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def make_sub_df_2(train, merge_df):\n    # Need to make week int string for join\n    train['Weeks'] = train['Weeks'].apply(str)\n    merge_df['Weeks'] = merge_df['Weeks'].apply(str)\n\n    new_sub = pd.merge(left=train, right=merge_df, on=['Patient','Weeks'], how='inner')\n\n    new_sub['Patient_Week'] = new_sub['Patient'].astype(str)+'_'+new_sub['Weeks'].astype(str)\n\n    new_sub = new_sub[['Patient_Week', 'FVC', 'Percent']]\n    new_sub = new_sub.rename(columns={'Percent': 'Confidence_actual', 'FVC': 'FVC_actual'})\n    \n    # change week back to int\n    train['Weeks'] = train['Weeks'].apply(int)\n    merge_df['Weeks'] = merge_df['Weeks'].apply(int)\n\n    return new_sub","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Final female and male sub dfs that are the second param to effnet_iter\n\nfemale_sub = make_sub_df_2(train, tempy_female)\nmale_sub = make_sub_df_2(train, tempy_male)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def effnet_iter(models, sub, train):\n    # sub needs to be with the Patient_Week column (female_sub or male_sub)\n    # train needs to be unique male_df or female_df\n    subs = []\n    for model in models:\n        q = 0.5\n#         sub = temp_fem.iloc[:3].copy()\n#         train = female_df.iloc[:3].copy()\n\n        A_test, B_test, P_test,W, FVC= {}, {}, {},{},{} \n        STD, WEEK = {}, {} \n        for p in tqdm(train.Patient.unique()):\n            x = [] \n            tab = [] \n            ldir = os.listdir(f'../input/osic-pulmonary-fibrosis-progression/train/{p}/')\n            for i in ldir:\n                if int(i[:-4]) / len(ldir) < 1.1 and int(i[:-4]) / len(ldir) > -0.1:\n                    x.append(get_img(f'../input/osic-pulmonary-fibrosis-progression/train/{p}/{i}')) \n                    tab.append(get_tab(train.loc[train.Patient == p, :])) \n            if len(x) <= 1:\n                continue\n            # tab is a list containing the csv data for all the images in a users folder\n            tab = np.array(tab) \n            # x is a list containing the img data for all the images in a users folder\n            x = np.expand_dims(x, axis=-1) \n            _a = model.predict([x, tab])\n            # a is the median of the predicted output _a\n            a = np.quantile(_a, q)\n\n            # We keep the median from the output of the NN\n            A_test[p] = a\n            # For the given patient, we take their FVC value then from that \n            # subtract the median times the week collected for that patient\n            B_test[p] = train.FVC.values[train.Patient == p] - a*train.Weeks.values[train.Patient == p]\n            # Keep the percent value found in the test df for the given patient\n            P_test[p] = train.Percent.values[train.Patient == p] \n            # Keep week info for given patient\n            WEEK[p] = train.Weeks.values[train.Patient == p]\n\n\n        # SO if we want to do it on a subset of data, \n        # we need a df that matches sub = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/sample_submission.csv')\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        _sub = sub.copy()\n        subs.append(_sub)\n    return subs","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fem_output = effnet_iter(models, female_sub.copy(), female_df.copy())","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fem_output[0]","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"male_output = effnet_iter(models, male_sub.copy(), male_df.copy())","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"male_output[0]","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Eval Metric\n\nThe error is thresholded at 1000 ml to avoid large errors adversely penalizing results, while the confidence values are clipped at 70 ml to reflect the approximate measurement uncertainty in FVC. The final score is calculated by averaging the metric across all test set Patient_Weeks threeperpatient. Note that metric values will be negative and higher is better.","metadata":{}},{"cell_type":"code","source":"def eval_metric(df):\n\n    sigma_clipped = [np.min([x, 70]) for x in df['Confidence_actual']]\n    temp = df['FVC_actual']-df['FVC']\n    delta = [np.min([np.abs(x), 1000]) for x in temp]\n\n    metric = []\n    for i in range(len(sigma_clipped)):\n        metric.append(((-1*np.sqrt(2)*delta[i])/sigma_clipped[i])-np.log(np.sqrt(2)*sigma_clipped[i]))\n        \n    return np.mean(metric)\n\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"female_score = eval_metric(fem_output[0])\nmale_score = eval_metric(male_output[0])","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"female_score, male_score","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"## evaluation metric function\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\n\n## default benchmark\n# male_ll_default = laplace_log_likelihood(male_output[0].FVC_actual, np.mean(male_output[0].FVC_actual), np.std(male_output[0].FVC_actual))\n# female_ll_default = laplace_log_likelihood(fem_output[0].FVC_actual, np.mean(fem_output[0].FVC_actual), np.std(fem_output[0].FVC_actual))\n\nmale_ll_default = laplace_log_likelihood(male_output[0].FVC_actual, np.mean(male_output[0].FVC_actual), male_output[0].Confidence_actual)\nfemale_ll_default = laplace_log_likelihood(fem_output[0].FVC_actual, np.mean(fem_output[0].FVC_actual), fem_output[0].Confidence_actual)\n\nmale_ll_pred = laplace_log_likelihood(male_output[0].FVC_actual, male_output[0].FVC, male_output[0].Confidence)\nfemale_ll_pred = laplace_log_likelihood(fem_output[0].FVC_actual, fem_output[0].FVC, fem_output[0].Confidence)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"female_ll_default, male_ll_default","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"female_ll_pred, male_ll_pred","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Original Implementation","metadata":{}},{"cell_type":"code","source":"subs = []\nfor model in models:\n\n    q = 0.5\n\n    sub = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/sample_submission.csv') \n    test = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/test.csv') \n    A_test, B_test, P_test,W, FVC= {}, {}, {},{},{} \n    STD, WEEK = {}, {} \n    for p in tqdm(test.Patient.unique()):\n        x = [] \n        tab = [] \n        ldir = os.listdir(f'../input/osic-pulmonary-fibrosis-progression/test/{p}/')\n        for i in ldir:\n            if int(i[:-4]) / len(ldir) < 1.1 and int(i[:-4]) / len(ldir) > -0.1:\n                x.append(get_img(f'../input/osic-pulmonary-fibrosis-progression/test/{p}/{i}')) \n                tab.append(get_tab(test.loc[test.Patient == p, :])) \n        if len(x) <= 1:\n            continue\n        # tab is a list containing the csv data for all the images in a users folder\n        tab = np.array(tab) \n        # x is a list containing the img data for all the images in a users folder\n        x = np.expand_dims(x, axis=-1) \n        _a = model.predict([x, tab])\n        # a is the median of the predicted output _a\n        a = np.quantile(_a, q)\n\n        # We keep the median from the output of the NN\n        A_test[p] = a\n        # For the given patient, we take their FVC value then from that \n        # subtract the median times the week collected for that patient\n        B_test[p] = test.FVC.values[test.Patient == p] - a*test.Weeks.values[test.Patient == p]\n        # Keep the percent value found in the test df for the given patient\n        P_test[p] = test.Percent.values[test.Patient == p] \n        # Keep week info for given patient\n        WEEK[p] = test.Weeks.values[test.Patient == p]\n        \n\n    # SO if we want to do it on a subset of data \n    # we need a df that matches sub = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/sample_submission.csv') \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)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Averaging Predictions","metadata":{}},{"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[0][\"FVC\"] * (1/N)\n    sub[\"Confidence\"] += subs[0][\"Confidence\"] * (1/N)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sub.head()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sub[[\"Patient_Week\",\"FVC\",\"Confidence\"]].to_csv(\"submission_img.csv\", index=False)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"img_sub = sub[[\"Patient_Week\",\"FVC\",\"Confidence\"]].copy()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Osic-Multiple-Quantile-Regression","metadata":{}},{"cell_type":"code","source":"ROOT = \"../input/osic-pulmonary-fibrosis-progression\"\nBATCH_SIZE=128\n\ntr = pd.read_csv(f\"{ROOT}/train.csv\")\ntr.drop_duplicates(keep=False, inplace=True, subset=['Patient','Weeks'])\nchunk = pd.read_csv(f\"{ROOT}/test.csv\")\n\nprint(\"add infos\")\nsub = pd.read_csv(f\"{ROOT}/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\")","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tr['WHERE'] = 'train'\nchunk['WHERE'] = 'val'\nsub['WHERE'] = 'test'\ndata = tr.append([chunk, sub])","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(tr.shape, chunk.shape, sub.shape, data.shape)\nprint(tr.Patient.nunique(), chunk.Patient.nunique(), sub.Patient.nunique(), \n      data.Patient.nunique())\n#","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"data['min_week'] = data['Weeks']\ndata.loc[data.WHERE=='test','min_week'] = np.nan\ndata['min_week'] = data.groupby('Patient')['min_week'].transform('min')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"base = data.loc[data.Weeks == data.min_week]\nbase = base[['Patient','FVC']].copy()\nbase.columns = ['Patient','min_FVC']\nbase['nb'] = 1\nbase['nb'] = base.groupby('Patient')['nb'].transform('cumsum')\nbase = base[base.nb==1]\nbase.drop('nb', axis=1, inplace=True)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"data = data.merge(base, on='Patient', how='left')\ndata['base_week'] = data['Weeks'] - data['min_week']\ndel base","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"COLS = ['Sex','SmokingStatus'] #,'Age'\nFE = []\nfor col in COLS:\n    for mod in data[col].unique():\n        FE.append(mod)\n        data[mod] = (data[col] == mod).astype(int)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#\ndata['age'] = (data['Age'] - data['Age'].min() ) / ( data['Age'].max() - data['Age'].min() )\ndata['BASE'] = (data['min_FVC'] - data['min_FVC'].min() ) / ( data['min_FVC'].max() - data['min_FVC'].min() )\ndata['week'] = (data['base_week'] - data['base_week'].min() ) / ( data['base_week'].max() - data['base_week'].min() )\n#data['percent'] = (data['Percent'] - data['Percent'].min() ) / ( data['Percent'].max() - data['Percent'].min() )\nFE += ['age','week','BASE']","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tr = data.loc[data.WHERE=='train']\nchunk = data.loc[data.WHERE=='val']\nsub = data.loc[data.WHERE=='test']\ndel data","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tr.shape, chunk.shape, sub.shape","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"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 make_model(nh):\n    z = L.Input((nh,), name=\"Patient\")\n    x = L.Dense(100, activation=\"relu\", name=\"d1\")(z)\n    x = L.Dense(100, activation=\"relu\", name=\"d2\")(x)\n    #x = L.Dense(100, activation=\"relu\", name=\"d3\")(x)\n    p1 = L.Dense(3, activation=\"linear\", name=\"p1\")(x)\n    p2 = L.Dense(3, activation=\"relu\", name=\"p2\")(x)\n    preds = L.Lambda(lambda x: x[0] + tf.cumsum(x[1], axis=1), \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.8), 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","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"y = tr['FVC'].values\nz = tr[FE].values\nze = sub[FE].values\nnh = z.shape[1]\npe = np.zeros((ze.shape[0], 3))\npred = np.zeros((z.shape[0], 3))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"net = make_model(nh)\nprint(net.summary())\nprint(net.count_params())","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"NFOLD = 2 # originally 5\nkf = KFold(n_splits=NFOLD)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\ncnt = 0\nEPOCHS = 600\nfor tr_idx, val_idx in kf.split(z):\n    cnt += 1\n    print(f\"FOLD {cnt}\")\n    net = make_model(nh)\n    net.fit(z[tr_idx], y[tr_idx], batch_size=BATCH_SIZE, epochs=EPOCHS, \n            validation_data=(z[val_idx], y[val_idx]), verbose=0) #\n    print(\"train\", net.evaluate(z[tr_idx], y[tr_idx], verbose=0, batch_size=BATCH_SIZE))\n    print(\"val\", net.evaluate(z[val_idx], y[val_idx], verbose=0, batch_size=BATCH_SIZE))\n    print(\"predict val...\")\n    pred[val_idx] = net.predict(z[val_idx], batch_size=BATCH_SIZE, verbose=0)\n    print(\"predict test...\")\n    pe += net.predict(ze, batch_size=BATCH_SIZE, verbose=0) / NFOLD","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sigma_opt = mean_absolute_error(y, pred[:, 1])\nunc = pred[:,2] - pred[:, 0]\nsigma_mean = np.mean(unc)\nprint(sigma_opt, sigma_mean)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"idxs = np.random.randint(0, y.shape[0], 100)\nplt.plot(y[idxs], label=\"ground truth\")\nplt.plot(pred[idxs, 0], label=\"q25\")\nplt.plot(pred[idxs, 1], label=\"q50\")\nplt.plot(pred[idxs, 2], label=\"q75\")\nplt.legend(loc=\"best\")\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(unc.min(), unc.mean(), unc.max(), (unc>=0).mean())","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.hist(unc)\nplt.title(\"uncertainty in prediction\")\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sub.head()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# PREDICTION\nsub['FVC1'] = 1.*pe[:, 1]\nsub['Confidence1'] = pe[:, 2] - pe[:, 0]\nsubm = sub[['Patient_Week','FVC','Confidence','FVC1','Confidence1']].copy()\nsubm.loc[~subm.FVC1.isnull()].head(10)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"subm.loc[~subm.FVC1.isnull(),'FVC'] = subm.loc[~subm.FVC1.isnull(),'FVC1']\nsigma_mean = 60\nif sigma_mean<sigma_mean:\n    subm['Confidence'] = sigma_opt\nelse:\n    subm.loc[~subm.FVC1.isnull(),'Confidence'] = subm.loc[~subm.FVC1.isnull(),'Confidence1']","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"subm.head()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"subm.describe().T","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"otest = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/test.csv')\nfor i in range(len(otest)):\n    subm.loc[subm['Patient_Week']==otest.Patient[i]+'_'+str(otest.Weeks[i]), 'FVC'] = otest.FVC[i]\n    subm.loc[subm['Patient_Week']==otest.Patient[i]+'_'+str(otest.Weeks[i]), 'Confidence'] = 0.1","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"subm[[\"Patient_Week\",\"FVC\",\"Confidence\"]].to_csv(\"submission_regression.csv\", index=False)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"reg_sub = subm[[\"Patient_Week\",\"FVC\",\"Confidence\"]].copy()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Ensemble (Simple Blend)","metadata":{}},{"cell_type":"code","source":"img_sub","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"reg_sub","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"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)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df = df1[['Patient_Week']].copy()\ndf['FVC'] = (0.45*df1['FVC'] + 0.55*df2['FVC'])\ndf['Confidence'] = (0.45*df1['Confidence'] + 0.55*df2['Confidence'])\ndf.head()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df.to_csv('submission.csv', index=False)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Test Code JIC","metadata":{}},{"cell_type":"code","source":"tempy_hehe = pd.DataFrame(columns=['Patient', 'Weeks'])\n\nfor i in tqdm(female_df.Patient):\n    for j in range(-12,133+1):\n        temp_name = i + '_' + str(j)\n#         print(i + '_' + str(j))\n        tempy_hehe = tempy_hehe.append({'Patient':i, 'Weeks' : j}, ignore_index=True)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tempy_hehe = tempy_hehe.sort_values(by=['Weeks'], ascending=[True])\n\ntempy_hehe['Patient_Week'] = tempy_hehe['Patient'].astype(str)+'_'+tempy_hehe['Weeks'].astype(str)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tempy_hehe['FVC'] = 0 \ntempy_hehe['FVC'] = 0 ","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tempy_hehe[['Patient_Week']]","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# This was the original way I was getting the Sub df, \n# but its not good bc it doesnt get all the img files for a given patient\n# need to go to the directory for each patient, get all their img\n\ntemp_fem = pd.DataFrame(columns=['Patient_Week', 'FVC', 'Confidence', 'FVC_actual', 'Percent_actual'])\n\ntemp_fem['Patient_Week'] = female_df['Patient'].astype(str)+'_'+female_df['Weeks'].astype(str)\ntemp_fem['FVC_actual'] = female_df['FVC']\ntemp_fem['Percent_actual'] = female_df['Percent']\n\ntemp_male = pd.DataFrame(columns=['Patient_Week', 'FVC', 'Confidence', 'FVC_actual', 'Percent_actual'])\n\ntemp_male['Patient_Week'] = male_df['Patient'].astype(str)+'_'+male_df['Weeks'].astype(str)\ntemp_male['FVC_actual'] = male_df['FVC']\ntemp_male['Percent_actual'] = male_df['Percent']\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}