{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"markdown","source":"## Objective:\nThe aim of this notebook is to document the experiments that I did during this competition. I had a ton of fun doing it. Yet it wasn't always fun & there were times when nothing I did seemed to improve the Public LB and I would be so frustrated. The metric was puzzling enough for me so I decided that I am not going to venture into the world of deep learning until I thoroughly completed my exploration with ML. However even this seemed a difficult task, so I decided to stick simply with `LinearRegression`.\n\nChoosing such a simple model simplified things greatly for me & provided me several advantanges:: \n1. There was no overfitting to the dataset. \n2. No parameters that needed fine tuning.\n3. It forces you to hand craft clever features that have more predictive potential.\n4. It is a white box algorithm, meaning that you can literally see how a model is making its predictions.\n\nUnfortunately the competition came to an end just when I was about to start using the scan data. Or maybe it was the other way round, I frantically tried incorporating and extracting features from the scan data. I have however documented my ideas as to how I might have used the scan data. Let me know if those ideas would indeed have worked.\n\n**It is my earnest hope that you derive the same joy that I experienced during my explorations by reading this Notebook. Happy reading!**\n\n## Import necessary modules, load the data:\n"},{"metadata":{},"cell_type":"markdown","source":"The snippet below is to install the GDCM module. Certain scans (Two of them) couldnt be opened simply via the `pydicom` module. We are unsure if any of the test scans contained similar scans. So it's better that we install them to be on the safer side."},{"metadata":{"trusted":true},"cell_type":"code","source":"!dpkg -i ../input/python3gdcm/build_1-1_amd64.deb\n!apt-get install -f\n\n!cp /usr/local/lib/gdcm.py /opt/conda/lib/python3.7/site-packages/.\n!cp /usr/local/lib/gdcmswig.py /opt/conda/lib/python3.7/site-packages/.\n!cp /usr/local/lib/_gdcmswig.so /opt/conda/lib/python3.7/site-packages/.\n!cp /usr/local/lib/libgdcm* /opt/conda/lib/python3.7/site-packages/.\n!ldconfig","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Import modules we would be frequently using. We would import many more along the way as they are needed."},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport tensorflow as tf\nimport gdcm\nimport pydicom\nimport glob\n\nfrom tqdm import tqdm\n%matplotlib inline","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# kaggle\nmain_dir = \"../input/osic-pulmonary-fibrosis-progression\"\n\n!ls {main_dir}","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_files = tf.io.gfile.glob(main_dir+\"/train/*/*\")\ntest_files = tf.io.gfile.glob(main_dir+\"/test/*/*\")\nsample_sub = pd.read_csv(main_dir + \"/sample_submission.csv\")\ntrain = pd.read_csv(main_dir + \"/train.csv\")\ntest = pd.read_csv(main_dir + \"/test.csv\")\n\nprint (\"Number of train patients: {}\\nNumber of test patients: {:4}\"\n       .format(train.Patient.nunique(), test.Patient.nunique()))\n\nprint (\"\\nTotal number of Train patient records: {}\\nTotal number of Test patient records: {:6}\"\n       .format(len(train_files), len(test_files)))\n\ntrain.shape, test.shape, sample_sub.shape","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Let's aquaint ourselves with reading, extracting DCM files using pydicom:"},{"metadata":{"trusted":true},"cell_type":"code","source":"# pydicom files are read with dcmread\ntemp = pydicom.dcmread(train_files[0])\ntype(temp)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# displaying a few items from the dicom file\n# elements are stored as tuple key pairs\nprint ('\\n'.join(str(temp).split(\"\\n\")[:15]))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Str function on `pydicom.dataset.FileDataset` returns a string representation.\nThe left part shows the keys (TAG NUMBERS or DICOM KEYWORDS) right part display the coresponding values."},{"metadata":{"trusted":true},"cell_type":"code","source":"# tuple keys can be obtained using .keys()\n# since `pydicom.dataset.FileDataset` is a \n# wrapper for dictionary files\nlist(temp.keys())[:5]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# returns all NON PRIVATE elements\n# values can be accessed both ways\nprint (temp.dir()[:5])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# a value can be accessed using all these methods\n(\n    temp.BitsAllocated, \n    temp.get('BitsAllocated'),\n    temp.data_element(\"BitsAllocated\").value, \n    temp[(0x28, 0x100)].value, \n    temp.get([0x28, 0x100]).value\n)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# We need to use .value when accessing using \n# tuple keys.But what does it exactly return?\nkey = (0x08, 0x08)\nprint (\"Accessing by a tuple key returns a\", type(temp[key]).__name__)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# some common methods of all DataElement(s)\nprint (list(filter(lambda x: \"__\" not in x, dir(pydicom.DataElement))))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# VR -> Value representation. The Value Representation of a Data Element \n# describes the data type and format of that Data Element's Value\ntemp[key].VR","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# other possible Value representations in this file alone\nnp.unique(list(map(lambda x: x.VR, temp.iterall())))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Common VR's observed:\n1. CS -> Code string. \n2. DS -> Decimal String\n3. IS -> Integer String\n4. LO -> Long String\n5. OW -> Other word \n6. PN -> Person Name\n7. SH -> Short string\n8. UI -> Unique identifier\n9. US -> Unsigned short\n\nVisit [here](http://dicom.nema.org/medical/dicom/current/output/chtml/part05/sect_6.2.html) for an exhaustive list of all possible value representations."},{"metadata":{"trusted":true},"cell_type":"code","source":"# usually the neatly formated\n# version of the DICOM Keyword\ntemp[key].description()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# number of values in the dataElement can be accessed using .VM\nprint (\"Value: {}\\nContains {} elements\".format(temp['Modality'].value, temp['Modality'].VM))\nprint (\"\\nValue: {}\\nContains {} elements\".format(temp['ImageType'].value, temp['ImageType'].VM))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# is the data element empty\ntemp['PatientSex'].is_empty, temp['ImageType'].is_empty","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# access DICOM keyword from tag number\ntemp[(0x18, 0x1151)].keyword","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# .repval returns value as a string\nprint (\"{} is saved as {}\".format(temp['XRayTubeCurrent'].repval, type(temp['XRayTubeCurrent'].repval)))\nprint (\"{} is saved as {}\".format(temp['XRayTubeCurrent'].value, type(temp['XRayTubeCurrent'].value)))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# all Dicom files contain a meta header that\n# can be accessed as follows:\ntemp.file_meta","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# to obtain a subset of the dataset\n# belonging to a key, 0x28 would\n# return all items of key tuples \n# with first element as 0x28\ntemp.group_dataset(0x28)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# to obtain the image from the dicom file\n# use pixel array; Bone is more suited for\n# ct scan images\nplt.figure(figsize=(8, 8))\nplt.axis('off')\nplt.imshow(temp.pixel_array, cmap='bone');","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Let's begin exploring our dataset:"},{"metadata":{"trusted":true},"cell_type":"code","source":"train.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# test has just 5 rows\ntest","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# we have to predict FVC & CONFIDENCE values\n# for these patients from week -12 to 133\n# (145 weeks each patient)\nsample_sub.tail()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# any missing values?\ntrain.isna().sum().any(), test.isna().sum().any()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# checking data type of each col\ntrain.info(null_counts=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"I suspect that some values in the train dataset are redundant, such as `Age`, `Sex` & `SmokingStatus`. Let's verify if we are correct:"},{"metadata":{"trusted":true},"cell_type":"code","source":"# columns that don't change with time for each patient\n(train.groupby('Patient').nunique() != 1).sum() == 0","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We are indeed correct. The above mentioned columns don't change over time at all, despite a very long time duration. Age of a person, Smoking Status & Sex (obviously) remains unchanged over time.\n\nAlthough we would expect the Age of a person to vary over a span of 145 weeks, this doesn't seem to happen.\n\n#### How many unique values have we got for each column?"},{"metadata":{"trusted":true},"cell_type":"code","source":"train.nunique()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"ages = train.groupby('Patient').Age.head(1)\nprint (\"Max Patient Age: {}\\nMin Patient Age: {}\".format(ages.min(), ages.max()))\nax = ages.plot(kind='hist', bins=50, edgecolor='red', color='y', figsize=(15, 5), xticks=range(49, 89))\nages.plot(kind='kde', ax=ax, xlim=(47, 90), color='w', secondary_y=True);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Observations:\n- Majority of the patients lie between the Age group of 60 - 75. \n- Has a guassian distribution."},{"metadata":{"trusted":true},"cell_type":"code","source":"f, ax = plt.subplots(figsize=(15, 5), ncols=2)\n\ntrain.groupby('Patient').SmokingStatus.head(1).value_counts().plot(\n    kind='pie', ax=ax[0], autopct=lambda x: str(int(x))+\"%\", \n    title='Smoking Status Pie chart', \n    colors=['orange', 'blue', 'green'])\n\ntrain.groupby('Patient').Sex.head(1).value_counts().plot(\n    kind='pie', ax=ax[1], autopct=lambda x: str(int(x))+\"%\", \n    title='Sex pie chart', colors=['red', 'blue']);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Observations:\n- A majority of the patients are ex smokers. \n- Most of the patients are male (~ 80%)\n\nHow does gender vary across smoking status?"},{"metadata":{"trusted":true},"cell_type":"code","source":"train.groupby(['SmokingStatus', 'Sex'])['Patient'].nunique().unstack().plot(\n    kind='bar', stacked=True, figsize=(10, 6), yticks=range(0, 130, 10),\n    rot=0, title='Gender Across Smoking Status');","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Observations:\n- Most of the Female patients have never smoked before\n- Most of the Male patients are ex smokers. \n- Very few people belong to the category of 'CurrentlySmoking'.\n\nLet's now see how the other features vary across the different genders:"},{"metadata":{"trusted":true},"cell_type":"code","source":"train.groupby(\"Sex\").agg(['min', 'max', 'mean', 'std']).drop(\"Age\", axis=1)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Some insights:\n- Males biologically have much higher FVC on average. \n- Even though mean FVC for males is greater, percent mean is actually lower for males than females. \n- Females are expected to be healthier? Females admitted are relatively healthier to males?\n\n#### How other features change with respect to age:"},{"metadata":{"trusted":true},"cell_type":"code","source":"f, ax = plt.subplots(nrows=2, figsize=(15, 10))\n\nsc = ax[0].scatter(\n    'Age', 'FVC', c=train.Sex.map({'Male': 0, 'Female': 1}), \n    s=(train.Weeks+5), data=train, cmap='brg_r', alpha=0.5)\n\nax[0].set(\n    xlabel='Age', ylabel='FVC', xticks=range(48, 90), \n    title=\"Age Vs FVC\")\nax[0].legend(sc.legend_elements()[0], ['Male', 'Female'])\n\nsc = ax[1].scatter(\n    'Age', 'Percent', c=train.Sex.map({'Male': 0, 'Female': 1}), \n    s=(train.Weeks+5), data=train, cmap='brg_r', alpha=0.5)\n\nax[1].set(\n    xlabel='Age', ylabel='FVC', xticks=range(48, 90), \n    title=\"Age Vs Percent\");\nax[1].legend(sc.legend_elements()[0], ['Male', 'Female'])\n\nf.suptitle(\"Across Different Genders\")\nf.tight_layout(rect=[0, 0.03, 1, 0.95]);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Observations:\n1. For any age group, FVC for males is higher than FVC for females as we had already observed\n2. Whereas, this distinction is not clear in case of Percent.\n\nLet's see statistically how the continous features vary across the categorical features:"},{"metadata":{"trusted":true},"cell_type":"code","source":"(train\n .groupby([\"SmokingStatus\"])[['Weeks', 'FVC', 'Percent']].agg({\n     \"Weeks\": \"count\",\n     \"FVC\": ['min', 'mean', 'max'], \n     \"Percent\": ['min', 'mean', 'max']})\n .rename({\"Weeks\": \"Cumulative Records\"}, axis=1))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Strangely enough, the current smokers have greater FVC & percent than other categories, as opposed to what we would expect. \n\nPossible explaination: Those who have come under our data might be there due a naturally poor FVC or be there *because* of smoking. Put in other words, *those who have never smoked and those who have naturally high FVC wouldn't be having much respiratory ailments or require diagnosis.*\n\nLet's now divide the above DataFrame by adding another category: Sex."},{"metadata":{"trusted":true},"cell_type":"code","source":"(train\n .groupby([\"SmokingStatus\", 'Sex'])[['Weeks', 'FVC', 'Percent']].agg({\n     \"Weeks\": \"count\",\n     \"FVC\": ['min', 'mean', 'max'], \n     \"Percent\": ['min', 'mean', 'max']})\n .rename({\"Weeks\": \"Cumulative Records\"}, axis=1))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"* Now the mean FVC's of those who have never smoked and the current smokers are a little closer. Although the current smokers still have higher FVC's, the difference isn't atleast as great as before. Maybe Smoking indeed increases FVC.\n\n* The percentage feature however is still unchanged.\n\n* Females who currently smoke have a much high Percent expectency than males. \n\n* The range of FVC for current female smokers are also within a much narrow window of 2700 to 2975. This is the case since, there are only two female patients who are current smokers:"},{"metadata":{"trusted":true},"cell_type":"code","source":"# females who smoke\ntrain.loc[(train.Sex == 'Female') & (train.SmokingStatus == 'Currently smokes'), 'Patient'].nunique()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"General FVC trend we observe: <br>\n1. FVC's for males: Never Smoked < Ex-Smokers < Current Smokers<br>\n2. FVC's for females: Never Smoked < Ex-smokers < Current Smokers.\n\n#### Let's calculate see the progression of FVC over the course of treatment for various categories:"},{"metadata":{"trusted":true},"cell_type":"code","source":"# we use savgol filter to smooth curves\nfrom scipy.signal import savgol_filter\n\ndef display_FVC_progress(data, title, smooth=True, drop=1, median=True):\n    \n    agg = ['count', 'min', 'median', 'max']\n    if not median:\n        agg.remove(\"median\")\n\n    temp = data.groupby('Weeks')[['FVC']].agg(agg)\n    temp = temp[temp['FVC']['count'] > drop].drop((\"FVC\", 'count'), axis=1)\n\n    if smooth:\n        # smoothing out the curves\n        temp['FVC', 'max'] = savgol_filter(temp['FVC', 'max'], 9, 3)\n        temp['FVC', 'min'] = savgol_filter(temp['FVC', 'min'], 9, 3)\n\n    ax = temp.plot(\n        figsize=(15, 5), \n        title=f'Variation & progress of FVC over the Weeks ({title})', \n        legend=True, xticks=range(-10, 150, 5));\n\n    ax.fill_between(temp.index, temp['FVC', 'max'], temp['FVC', 'min'], color='green');","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# FVC Progress across all patients\ndisplay_FVC_progress(train, 'All Categories')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"- We can observe that intially the variation is high but over the weeks, variation between max and min FVC keeps decreasing. Confidence value could initially be low and later on increased.\n\n- The median values are usually closer to min than max (again indicates the skew). *It's safer to predict values closer to min than to max.*\n\n- There is a slight declining trend of FVC"},{"metadata":{"trusted":true},"cell_type":"code","source":"# FVC progress across Sex\ndisplay_FVC_progress(train.loc[train.Sex == 'Male'], 'Only Males', drop=1, smooth=True, median=False)\ndisplay_FVC_progress(train.loc[train.Sex == 'Female'], 'Only Females', drop=1, smooth=True, median=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"- The trend for males is slightly declining. While it is slightly improving for females.\n- The range appears more chaotic for females, probably due to low number of female patient data available."},{"metadata":{"trusted":true},"cell_type":"code","source":"# FVC progress across SmokingStatus\ndisplay_FVC_progress(train.loc[train.SmokingStatus == 'Ex-smoker'], 'Cat: Ex Smokers', median=False)\ndisplay_FVC_progress(train.loc[train.SmokingStatus == 'Currently smokes'], 'Cat: Current Smokers', median=False)\ndisplay_FVC_progress(train.loc[train.SmokingStatus == 'Never smoked'], 'Cat: Never Smoked', median=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Ex smoker's FVC initially decreases and then flattens out. Whereas the trend for Never smoked and currently smoking is increasing. Although we need to remind ourselvs that there are relatively little Current smoker's  data available.\n\nNow that we have seen how the trend is for different groups, set's see how FVC progression of random patients from each group look like. But before we do that, let's answer a more simpler question: how many observations do we have for each patient?"},{"metadata":{"trusted":true},"cell_type":"code","source":"# min, max and mean patient observations\ntrain.groupby('Patient')['Weeks'].count().agg(['min', 'max', 'mean'])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"It varies from a minimum of 6 weeks to a maximum of 10 week observations, using which we have to predict for 145 weeks. "},{"metadata":{"trusted":true},"cell_type":"code","source":"# any repeating observations for the same week?\ntemp = train.groupby(\"Patient\")['Weeks'].agg(['count', 'nunique'])\ntemp = temp[temp['count'] != temp['nunique']].index\nlen(temp)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We have some patients for whom multiple FVC were taken for the same week, let's probe further:"},{"metadata":{"trusted":true},"cell_type":"code","source":"train[train.Patient == np.random.choice(temp)]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"All these patients with multiple observations taken for the same week occur either on the first week or on the last week of their clinical observation. \n\nMaybe they were taken for confirmation. Let's keep only that FVC observation that is closest to second Week's Observation:"},{"metadata":{"trusted":true},"cell_type":"code","source":"multiple = train.groupby(['Patient', \"Weeks\"])[['FVC', 'Percent']].agg(['min', 'max'], 1)\n\nmultiple = multiple[[('FVC', 'min'), (\"Percent\", 'min')]].where(\n    multiple['Percent'].diff().abs().shift(-1).apply(np.argmin, 1) == 0, \n    multiple[[('FVC', 'max'), (\"Percent\", 'max')]].values)\n\nmultiple.columns = multiple.columns.droplevel(1)\nmultiple = multiple.rename_axis(None, axis=1).reset_index()\n\ntrain = multiple.merge(\n    train[['Patient', 'Age', 'Sex', 'SmokingStatus']].groupby(\"Patient\").head(1),\n    on='Patient')\n\n# rain check\ntrain[train.Patient == np.random.choice(temp)]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### FVC progression samples of patients from various groups:"},{"metadata":{"trusted":true},"cell_type":"code","source":"# to increase choices, set replace as false in np.random.choice\n# since there are very few female patient's available\nchoice = 2 # two each\n\ntemp = (train.groupby(['Sex', 'SmokingStatus'])['Patient']\n        .apply(lambda x: np.random.choice(np.unique(x), size=choice, replace=False))\n        .reset_index())\n\nf, ax = plt.subplots(ncols=choice, nrows=6, figsize=(20, 30))\nfor i, sex, status, patients in temp.itertuples():\n    \n    for j in range(choice):\n    \n        (train\n         .loc[train.Patient == patients[j], ['FVC', 'Weeks']]\n         .set_index('Weeks')\n         .plot(ax=ax[i][j], title=f\"{sex} Patient\\n{status}\", legend=False)\n        )\n        \nf.tight_layout();","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The metric used in this competition is `Laplace log likelihood`. It works as follows:\n\n1. Confidence values smaller than 70 are clipped. \n\n$σ_{clipped}=max(σ,70)$\n\n2. Errors greater than 1000 are also clipped in order to avoid large errors. \n\n$Δ=min(|FVC_{true}−FVC_{predicted}|,1000)$\n\n3. Finally the metric is defined as: \n\n$metric=−\\dfrac{\\sqrt2Δ}{σ_{clipped}}−ln(\\sqrt2σ_{clipped})$\n\nWe simply rewrite the code from [here](https://www.kaggle.com/gunesevitan/osic-pulmonary-fibrosis-progression-eda) in numpy to tensorflow. *I wrote in tensorflow so that it would be easy to reuse the code in tf models. However my plans changed and I never really used TF much.*"},{"metadata":{"trusted":true},"cell_type":"code","source":"def laplace_log_likelihood(y_true, y_pred, sigma=70):\n    # values smaller than 70 are clipped\n    sigma_clipped = tf.maximum(sigma, 70)\n\n    # errors greater than 1000 are clipped\n    delta_clipped = tf.minimum(tf.abs(y_true - y_pred), 1000)\n    \n    # type cast them suitably\n    delta_clipped = tf.cast(delta_clipped, dtype=tf.float32)\n    sigma_clipped = tf.cast(sigma_clipped, dtype=tf.float32)\n    \n    # score function\n    score = - tf.sqrt(2.0) * delta_clipped / sigma_clipped - tf.math.log(tf.sqrt(2.0) * sigma_clipped)\n    \n    return tf.reduce_mean(score)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# perfect score would be\nlaplace_log_likelihood(train['FVC'], train['FVC'], 70)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now that we have our scoring function ready, let's see how `confidence` affects the scores:"},{"metadata":{"trusted":true},"cell_type":"code","source":"# an impossibly high FVC (delta would be set to 1000)\n# varying confidence value to understand how it would affect the score\nhigh_delta = []\nzero_delta = []\nfor i in range(70, 2000):\n    high_delta.append(laplace_log_likelihood(train['FVC'], 9e5, sigma=i).numpy())\n    zero_delta.append(laplace_log_likelihood(train['FVC'], train['FVC'], sigma=i).numpy())\n    \n   \n# lets see how score varies with differing conf\nf, ax = plt.subplots(figsize=(20, 5), ncols=2)\nax[0].plot(high_delta)\nax[0].set(xlabel='Confidence', ylabel='Scores', title='$\\Delta = 1000$ (Incorrect Predictions)')\n\nax[1].plot(zero_delta)\nax[1].set(xlabel='Confidence', ylabel='Scores', title='$\\Delta = 0$ (Correct Predictions)');","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"In a hypothetical situation where $\\Delta$ value can only be 0 or 1000:\n\n1. Predict confidence value as 70 when you know you are correct\n2. Predict confidence value as 250 or greater when you know you are wrong.\n\n*The insight gained*: The closer we are to the actual FVC, stick with lower confidence values. The farther we from actual FVC, stick with large confidence values.\n\nLet's now create the most basic of all baselines, which would be to predict FVC such that all $\\Delta$ values are greater than 1000. Applying the knowledge we have with confidence let's set the confidence value as 250 for this baseline and see how it scores on the LB:"},{"metadata":{"trusted":true},"cell_type":"code","source":"sub = sample_sub.copy()\nsub.FVC = 9e3             # we are clear off the max FVC by atleast 1000\nsub.Confidence = 250      # Increasing Conf would directly increase the Score, Try yourself\nsub.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# sub.to_csv(\"confidence_check.csv\", index=False)\nsub.to_csv(\"conf_submission.csv\", index=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Submission Scores:\n- Public LB score: `-11.5248`\n- Private LB score: `-11.5248`\n\n\n#### Let's explore some less dumb, less naive ideas before we move on to more concrete ideas."},{"metadata":{"trusted":true},"cell_type":"code","source":"# binned ages for grouping\nages = pd.cut(train.Age, 10).cat.codes\n\ndumb_preds = [\n    (\"Sample Submission Idea\", 2000),\n    \n    (\"Min Scores\", train['FVC'].min()), \n    (\"25th Quantile Scores\", train['FVC'].quantile(0.25)), \n    (\"Median Scores\", train['FVC'].median()), \n    (\"75th Quantile Scores\", train['FVC'].quantile(0.75)), \n    (\"Max Scores\", train['FVC'].max()), \n    \n    (\"Mean Scores\", train['FVC'].mean()),\n    \n    (\"Weeks median\", train.groupby('Weeks')['FVC'].transform('median')), \n    (\"Binned Age median\", train.groupby([ages])['FVC'].transform('median')),\n    (\"SmokingStatus median\", train.groupby(['SmokingStatus'])['FVC'].transform('median')),\n    (\"Sex median\", train.groupby(['Sex'])['FVC'].transform('median')),\n    \n    (\"Sex-SmokingStatus median\", train.groupby(['Sex', 'SmokingStatus'])['FVC'].transform('median')),\n    (\"Age-Sex median\", train.groupby([ages, 'Sex'])['FVC'].transform('median')),\n    (\"Age-SmokingStatus median\", train.groupby([ages, 'SmokingStatus'])['FVC'].transform('median')),\n    (\"Weekly-Sex median\", train.groupby(['Weeks', 'Sex'])['FVC'].transform('median')),\n    (\"Weekly-Smoking median\", train.groupby(['Weeks', 'SmokingStatus'])['FVC'].transform('median')),\n    (\"Weekly-Age median\", train.groupby(['Weeks', ages])['FVC'].transform('median')),\n    \n    (\"Weekly-Sex-Smoking median\", train.groupby(['Weeks', 'Sex', 'SmokingStatus'])['FVC'].transform('median')),\n    (\"Weekly-Sex-Age median\", train.groupby([\"Weeks\", \"Sex\", ages])['FVC'].transform('median')),\n    (\"Weekly-Smoking-Age median\", train.groupby([\"Weeks\", \"SmokingStatus\", ages])['FVC'].transform('median')),\n]\n\n# fixed confidence\nsigma = 250\n\n# confidence increases slowly towards the end \n# Since we observed that range of FVC decreases \n# with increasing Weeks\nsigma_l = train.groupby('Patient')['Weeks'].transform(lambda x: np.linspace(225, 275, len(x))).values\n\n\nprint (\"Some Dumb Ideas & their Scores:\\n\")\nfor text, preds in dumb_preds:\n    \n    score = laplace_log_likelihood(train['FVC'], preds, sigma).numpy()\n    print (f\"\\t{text} with fixed conf {' ' * (29 - len(text))}: {score:-6.2f}\")\n    \n    score = laplace_log_likelihood(train['FVC'], preds, sigma_l).numpy()\n    print (f\"\\t{text} with conf swelling {' ' * (26 - len(text))}: {score:-6.2f}\\n\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"markdown","source":"Confidence swelling does very little to improve the score probably since we are using medians and means, it might work better when we predict with advanced models.\n\nWe also observe that weekly-Smoking-Age is able to obtain the best possible score. However, we need to ensure that every groupby in groupby has sufficient rows, otherwise median would return a much too biased score (Since we are directly obtaining our pred from FVCs). Let's ensure the same:"},{"metadata":{"trusted":true},"cell_type":"code","source":"# number of groups with more than 1 rows out of 1549 instances\n(train.groupby([\"Weeks\", \"SmokingStatus\", ages])['FVC'].count() != 1).sum()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's calculate for the other groups as well. We use mean to measure how reliable the score is:"},{"metadata":{"trusted":true},"cell_type":"code","source":"print (\"{} {}| {:^5} | {:^5} | {:^5} | {:^5}\\n{}\".format(\n    'Category (Patients Available)', \n    ' ' * (30 - len('Category (Patients Available)')), \n    'Min', 'Med', 'Max', \"Count\", \"=\"*65)\n)\n\nfor name, group in (\n    (\"Weekly\", train.groupby('Weeks')), \n    (\"Binned Age\", train.groupby([ages])),\n    (\"SmokingStatus\", train.groupby(['SmokingStatus'])),\n    (\"Sex\", train.groupby(['Sex'])),\n    (\"Sex-SmokingStatus\", train.groupby(['Sex', 'SmokingStatus'])),\n    (\"BAge-Sex\", train.groupby([ages, 'Sex'])),\n    (\"BAge-SmokingStatus\", train.groupby([ages, 'SmokingStatus'])),\n    (\"Weekly-Sex\", train.groupby(['Weeks', 'Sex'])),\n    (\"Weekly-SmokingStatus\", train.groupby(['Weeks', 'SmokingStatus'])),\n    (\"Weekly-BAge\", train.groupby(['Weeks', ages])),\n    (\"Weekly-Sex-Smoking\", train.groupby(['Weeks', 'Sex', 'SmokingStatus'])),\n    (\"Weekly-Sex-BAge\", train.groupby([\"Weeks\", \"Sex\", ages])),\n    (\"Weekly-Smoking-BAge\", train.groupby([\"Weeks\", \"SmokingStatus\", ages])),\n):\n    \n    mn, md, mx, cnt = group['Patient'].nunique().agg(['min', 'median', 'max', 'count']).tolist()\n    print (\"{} {}| {:5.1f} | {:5.1f} | {:5.1f} | {:5.1f}\".format(name, ' ' * (30 - len(name)), mn, md, mx, cnt))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Age is really good, however there are very few instances available for meaningful predictions. Same applies for almost all the categories. We would be simply copying the stats over to the test set this way. Let's continue however and see the top 3 ideas that work well:\n1. Weekly-Sex-Smoking Median: -8.06 (med: 2)\n2. Weekly-Sex Median: -8.35 (med: 4) **\n3. Weekly-Smoking Median: -8.58 (med: 4)\n\nAll these well performing ideas require data available for weeks for which we are to predict. This is a problem in those cases where weeks such as -12 need to be predicted. Let's choose some ideas which work better without needing week's data:\n\n1. Sex-SmokingStatus median:  -8.48 (med: 17.5)\n3. BAge-Sex median: -8.48 (med: 5.5)\n2. Sex median: -8.60 (med: 88) **\n\nOur submission idea: \n\nWe use Weekly-Sex median for those availble and for those where week data is not available we use the Sex median. We need to create the submission dataframe first:"},{"metadata":{"trusted":true},"cell_type":"code","source":"sub = sample_sub.Patient_Week.str.extract(\"(ID\\w+)_(\\-?\\d+)\").rename({0: \"Patient\", 1: \"Weeks\"}, axis=1)\nsub['Weeks'] = sub['Weeks'].astype(int)\nsub = pd.merge(sub, test[['Patient', 'Sex', 'SmokingStatus']], on='Patient')\nsub.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's now use the dataframe to apply our logic to get the baseline predictions on test data:"},{"metadata":{"trusted":true},"cell_type":"code","source":"week_temp = train.groupby([\"Weeks\", 'Sex'])['FVC'].median()\nsex_temp = train.groupby(['Sex'])['FVC'].median()\n\nfor index, week, sex in sub.iloc[:, 1:3].itertuples():\n    if (week, sex) in week_temp:\n        # we assume we are more accurate here\n        sub.loc[index, 'FVC'] = week_temp[week, sex]\n        sub.loc[index, 'Confidence'] = sigma\n    else:\n        # we assume we are less accurate here, boost confidence\n        sub.loc[index, 'FVC'] = sex_temp[sex]\n        sub.loc[index, 'Confidence'] = sigma + 100\n        \nsub.sample(5)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# swelling confidence as progress in the weeks\nsub[\"Patient_Week\"] = sub.Patient + \"_\" + sub.Weeks.astype(str)\nsub.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# naive submission\nsub[['Patient_Week', 'FVC', 'Confidence']].to_csv(\"pd_submission.csv\", index=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Submission Scores:\n- Public LB score: `-9.0917`\n- Private LB score: `-8.7322`\n\n*In hindsight, I can't beleive that this naive idea actually did well in Private LB*\n\nLet's now use the data we have on the train dataset itself to create a better **model** and see if we can perform better. We use only the `Weeks`, `Age`, `Sex`, `SmokingStatus`:\n\n*Note*: We refrain from using `Percent` since that data is not available for the test data (as of now)."},{"metadata":{"trusted":true},"cell_type":"code","source":"x = train[['Weeks', 'Age', 'Sex', 'SmokingStatus', 'Patient']].copy()\ny = train['FVC'].copy()\n\n# save the stats for future use\nstats = x.describe().T\n\n# one hot encoding\nx = pd.get_dummies(x, columns=['Sex', 'SmokingStatus'], drop_first=True)\n\n# scaling the numeric features\nfor col in ['Weeks', 'Age']:\n    x[col] = (x[col] - stats.loc[col, 'min']) / (stats.loc[col, 'max'] - stats.loc[col, 'min'])\n\n# how does it look?\nx.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.model_selection import cross_val_score, GroupKFold\nfrom sklearn.linear_model import LinearRegression\nfrom sklearn.metrics import make_scorer\n\n# we use groupkfold to ensure there is no data leakage\ncv = GroupKFold(n_splits=7)\n\n# create a simple scorer function that can alter the value of sigma\n# further we return an Sklearn Scorer to be able to fit inside\n# sklearn's models and pipelines\ndef l1(s):\n    def scorer_func(x, y, sigma=s):\n        return laplace_log_likelihood(x, y, sigma=s).numpy()\n    \n    return make_scorer(scorer_func, greater_is_better=False)\n\ncross_val_score(LinearRegression(), x.drop(\"Patient\", 1), y, cv=cv, groups=x.Patient, scoring=l1(70))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"It doesn't do that well, scores similar to our naive predictions ;( \n\nLet's see if we can do better by feature engineering from existing columns.\n\n1. `FVC` is the value we are trying to predict. Future or past FVC values have to be closer to the base line FVC of the same person\n2. The problem with `Percent` is that it is available for all entry on the train dataset, however only first week's percent is available in case of the test dataset. Maybe we could create an auxillary model that predicts the percent feature. This model could better help the main model that does the job of predicting FVC. *We will do this at a later stage.*"},{"metadata":{"trusted":true},"cell_type":"code","source":"x = train.copy()\ny = train['FVC'].copy()\n\n# base features that can be useful for making predictions\nx['base_Week'] = x.groupby('Patient')['Weeks'].transform('min')\nx['Base_FVC'] = x.groupby('Patient')['FVC'].transform('first')\n\n# save the stats for future use\nstats = x.describe().T\n\n# one hot encoding\nx = pd.get_dummies(x, columns=['Sex', 'SmokingStatus'], drop_first=True)\n\n# numeric columns to a list\nnum_cols = [\n    'Weeks', 'Age', 'base_Week', 'Base_FVC'\n]\n\n# scaling the numeric features\nfor col in num_cols:\n    x[col] = (x[col] - stats.loc[col, 'min']) / (stats.loc[col, 'max'] - stats.loc[col, 'min'])\n    \n# print out how well our features would do\nprint (x.corr()['FVC'].abs().sort_values(ascending=False)[1:])\n\n# drop Percent for now\nx.drop(['FVC', 'Percent'], axis=1, inplace=True)\nx.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# how does it perform now?\ncross_val_score(\n    LinearRegression(), x.drop(\"Patient\", 1), y, \n    cv=cv, groups=x.Patient, \n    scoring=l1(70)\n).mean()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# fit on the train dataset\nlr = LinearRegression().fit(x.drop(\"Patient\", 1), y)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This model performs much better than previous model (atleast only on the folds).\n\nLet's make our predictions on test set and see how it scores:"},{"metadata":{"trusted":true},"cell_type":"code","source":"x = (sub.drop(['Confidence', 'Patient_Week'], 1)\n     .merge(test[['Patient', 'Weeks', 'FVC', 'Age']], on='Patient')\n     .rename({\"Weeks_y\": \"base_Week\", \"FVC_y\": \"Base_FVC\", \"Weeks_x\": \"Weeks\"}, axis=1)\n     .drop(['Patient', 'FVC_x'], axis=1))\n\n# one hot encoding, We set drop_first as \n# false to ensure the test is same as train\nx = pd.get_dummies(x, columns=['Sex', 'SmokingStatus'])\n\n# # scaling the numeric features\nfor col in ['Weeks', 'Age', 'base_Week', 'Base_FVC']:\n    x[col] = (x[col] - stats.loc[col, 'min']) / (stats.loc[col, 'max'] - stats.loc[col, 'min'])\n\nx = x[['Weeks', 'Age', 'base_Week', 'Base_FVC', 'Sex_Male',\n       'SmokingStatus_Ex-smoker', 'SmokingStatus_Never smoked']]\n\nx.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# LR submission\nsub['FVC'] = lr.predict(x)\nsub[['Patient_Week', 'FVC', 'Confidence']].to_csv(\"LR_submission_improved.csv\", index=False)\n\nsub.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Submission Scores:\n- Public LB score: `-6.9135`\n- Private LB score: `-6.8646`\n\nOne of the major drawback with our approach is that while the train dataset we only predict forwards, we also need to be able to predict backwards in time for our test data. So instead of setting the base week as the first week's record of the patient, we could do better if we set the base week as somewhere along the 25% of records we  have for that patient.\n\nLet's create a function for that:"},{"metadata":{"trusted":true},"cell_type":"code","source":"def base_shift(data, q=50):\n    '''\n    Create base_Week, Base_FVC and Base_Percent for train\n    based on the given percentile. q=0 sets the base_Week\n    as the first week in the dataset.\n    '''\n    \n    x = data.copy()\n    \n    temp = (x.groupby(\"Patient\")\n            .apply(lambda x: x.loc[int(\n                np.percentile(x['Weeks'].index, q=q)\n            ), [\"Weeks\", \"FVC\", \"Percent\"]]))\n\n    temp.rename(\n        {\"Weeks\": \"Base_Week\", \n         \"FVC\": \"Base_FVC\", \n         \"Percent\": \"Base_Percent\"}, \n        axis=1, inplace=True)\n\n    # merge it with train data\n    x = x.merge(temp, on='Patient')\n\n    # create week offsets\n    x['Week_Offset'] = x['Weeks'] - x['Base_Week']\n    \n    return x","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.preprocessing import OneHotEncoder\nonehcenc = OneHotEncoder()\n\nfolds = 7\ndata = {}\n\n# do a train/val split for evaluation of the strategy\ntotal_patients = train.Patient.unique()\nnp.random.shuffle(total_patients)\nval_len = len(total_patients) // folds\n\nfor i in range(folds):\n    \n    # do a train/val split for evaluation of the strategy   \n    val_patients = total_patients[(i)*val_len:(i+1)*val_len]\n    train_patients = np.setdiff1d(total_patients, val_patients)\n\n    # shift base by 25%\n    x = base_shift(train[train.Patient.isin(train_patients)], q=25)\n\n    # one hot encoding\n    x = x.merge(\n        pd.DataFrame(\n            onehcenc.fit_transform(x[['Sex', 'SmokingStatus']]).todense(),\n            columns=[*np.concatenate(onehcenc.categories_)]),\n        left_index=True, right_index=True\n    )\n\n    # binned FVC does better?\n    x['Bin_base_FVC'] = pd.cut(x['Base_FVC'], bins=range(0, 7501, 500)).cat.codes / 15\n\n    # saving stats for future\n    stats = x.describe().T\n\n    # lets scale the numeric columns (We scale it with max possibe values)\n    num_cols = ['Weeks', 'Week_Offset', 'Base_Week', 'Age', \n                'Base_FVC', 'Percent', 'Base_Percent']\n    \n    for col in num_cols:\n        x[col] = (x[col] - stats.loc[col, 'min']) / (stats.loc[col, 'max'] - stats.loc[col, 'min'])\n\n    to_drop = (\n        [\"FVC\", 'Percent', 'Sex', 'SmokingStatus']\n\n        + [\n    #         \"Base_FVC\", \n    #         'Base_Week', \n    #         'Weeks', \n    #         'Bin_base_FVC', \n    #         'Base_Percent'\n        ] \n    )\n\n    y = x['FVC'].dropna()\n    x = x.drop(to_drop, axis=1)\n    x_cols = x.columns\n    \n    x_val = base_shift(train[train.Patient.isin(val_patients)], q=5)\n\n    x_val = x_val.merge(\n        pd.DataFrame(\n            onehcenc.transform(x_val[['Sex', 'SmokingStatus']]).todense(),\n            columns=[*np.concatenate(onehcenc.categories_)]),\n        left_index=True, right_index=True\n    )\n\n    x_val['Bin_base_FVC'] = pd.cut(x_val['Base_FVC'], bins=range(0, 7501, 500)).cat.codes / 15\n\n    for col in num_cols:\n        if col in x_val.columns:\n            x_val[col] = (x_val[col] - stats.loc[col, \"min\"]) / (stats.loc[col, 'max'] - stats.loc[col, 'min'])\n\n    y_val = x_val['FVC'].dropna()        \n    x_val = x_val.drop(to_drop, axis=1, errors='ignore')\n    x_val = x_val[x_cols]\n    \n    # saving the folds\n    data[i] = (x, y, x_val, y_val)\n    \n\nprint (\"Total Number of Folds: \" + str(folds))\n\n# how does it look?\ndata[0][0].head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Some of the features on the training set seems rather useless. So let's create a simple pipeline to select only the most prominent of the features:"},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.feature_selection import SelectKBest, f_regression\nfrom sklearn.pipeline import Pipeline\nfrom sklearn.model_selection import GridSearchCV\n\ngrid_params = {\"SFromModel__k\": range(1, data[0][0].shape[1]+1)}\n\ntemp = Pipeline(\n    [(\"SFromModel\", SelectKBest(score_func=f_regression)),\n    (\"Model\", LinearRegression())])\n\nscores = []\n\nfor i in range(folds):\n    x, y, _, _ = data[i]\n    grid = GridSearchCV(temp, param_grid=grid_params, n_jobs=-1, cv=cv, scoring=l1(70))\n    grid.fit(x.drop(\"Patient\", 1), y, groups=x.Patient)\n\n    scores.append((grid.best_params_, grid.best_score_))\n    \nscores = sorted(scores, key=lambda x: x[1])\nbest_params = scores[0][0]\nbest_score = scores[0][1]\n\nprint (\"Mean Train Score: {:.2f}\\nBest Score: {:10.2f}\\nBest params: {:9}\".format(\n    sum(scores[i][1] for i in range(folds)) / folds,\n    best_score, best_params['SFromModel__k']))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# choosing the best parameter from above\nmodel = Pipeline(\n    [(\"SFromModel\", SelectKBest(score_func=f_regression, k=best_params['SFromModel__k'])),\n    (\"Model\", LinearRegression())])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's also see the value of sigma which scores the best with this model:"},{"metadata":{"trusted":true},"cell_type":"code","source":"# best score placeholder\nbest_score = (0, np.inf, np.inf)\n\nfor i in range(50, 1500, 50):\n    \n    temp = cross_val_score(model, x.drop(\"Patient\", 1), y, cv=cv, groups=x.Patient, scoring=l1(i))\n    if temp.mean() < best_score[1]:\n        best_score = i, temp.mean(), temp.std()\n        \nbest_score","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Before we get excited about this idea's score, let's verify with our validation dataset:"},{"metadata":{"trusted":true},"cell_type":"code","source":"scores = {}\n\nfor i in range(folds):\n    x, y, x_val, y_val = data[i]\n    \n    model.fit(x.drop(\"Patient\", 1), y)\n    scores[70] = scores.get(70, []) + [- laplace_log_likelihood(\n        y_val, model.predict(x_val.drop(\"Patient\", 1)), sigma=70)]\n    \n    scores[best_score[0]] = scores.get(best_score[0], []) + [- laplace_log_likelihood(\n        y_val, model.predict(x_val.drop(\"Patient\", 1)), sigma=best_score[0])]\n    \nprint (\"Evaluation Scores Mean: {:.2f} @  70 Confidence\\n\\\nEvaluation Scores Mean: {:.2f} @ {} Confidence\".format(\n        np.mean(scores[70]), \n        np.mean(scores[best_score[0]]), \n        best_score[0]\n    ))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Features that were selected can be visualized. This is pro of ML algorithms:"},{"metadata":{"trusted":true},"cell_type":"code","source":"# features that were chosen (ascending order of importance)\ntemp = x.columns.drop(\"Patient\").values[np.argsort(model.steps[0][1].scores_)].tolist()\nprint (\"Those that were picked :\", temp[-best_params['SFromModel__k']:])\nprint (\"Those that were dropped:\", temp[:best_params['SFromModel__k']])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Previously our model made a cross_val_score of *7.66*, this made a score lesser than that. This idea hasn't helped our model. Although it performs well in the train dataset, it doesn't seem to pick any of the week based feature and as such its no surprise that this model would perform poorly on the test set.\n\nLet's proceed with predictions and submit it, just since.. well simply because there's no real limit to making submissions at this point of the competition."},{"metadata":{"trusted":true},"cell_type":"code","source":"# merge the test dataset as well to be able to handle 1hC\nx = sub[['Patient', 'Weeks']].merge(\n    test.rename({\"Weeks\": \"Base_Week\", \n                 \"FVC\": \"Base_FVC\", \n                 \"Percent\": \"Base_Percent\"}, axis=1), \n    on='Patient')\n\n# create week offsets\nx['Week_Offset'] = x['Weeks'] - x['Base_Week']\n\n# one hot encoding\nx = x.merge(\n    pd.DataFrame(\n        onehcenc.transform(x[['Sex', 'SmokingStatus']]).todense(),\n        columns=[*np.concatenate(onehcenc.categories_)]),\n    left_index=True, right_index=True\n)\n\n# binned FVC does better?\nx['Bin_base_FVC'] = pd.cut(x['Base_FVC'], bins=range(0, 7501, 500)).cat.codes / 15\n\n# lets scale the numeric columns (We scale it with max possibe values)\nfor col in num_cols:\n    if col in x.columns:\n        x[col] = (x[col] - stats.loc[col, \"min\"]) / (stats.loc[col, 'max'] - stats.loc[col, 'min'])\n\nx = x.drop(to_drop, axis=1, errors='ignore')\n\nx = x[x_cols]\n\nx.drop(\"Patient\", 1).head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# LR submission\nsub['FVC'] = model.predict(x.drop(\"Patient\", 1))\nsub['Confidence'] = best_score[0]\n\n# final touches before submission: Simply copy paste those rows\n# we already know from the test dataset given with conf 70\nfor i in range(len(test)):\n    sub.loc[sub['Patient_Week']==test.Patient[i]+'_'+str(test.Weeks[i]), 'FVC'] = test.FVC[i]\n    sub.loc[sub['Patient_Week']==test.Patient[i]+'_'+str(test.Weeks[i]), 'Confidence'] = 70\n\n# submissions to sub file\nsub[['Patient_Week', 'FVC', 'Confidence']].to_csv(\"lr_submission_base_week_shifted.csv\", index=False)\n\n# how does it look?\nsub.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Submission Scores:\n- Public LB score: `-7.1386`\n- Private LB score: `-7.0275`\n\n*This idea has not helped our model. Also note the static FVC values xD.*\n\nOur same idea had been already deployed in [this](https://www.kaggle.com/yasufuminakama/osic-lgb-baseline) notebook. However the author of this kernel has suggested a much better idea to create multiple base_weeks instead of choosing one in the 25% or in the start, this way we would get around 7x more data to train our model with. \n\nIsn't the adage goes, \"Quality data over fancier algorithms\"? So, Let's create a function to do that:"},{"metadata":{"trusted":true},"cell_type":"code","source":"def multi_baseweek_frame(data, display=True):\n    '''Function to return multiple base week frames -> instead of creating one base week,\n    creates several to help our model learn better and predict past and future data better.\n    \n    Also this function works much faster than the original function from Y.Nakama\n    since it uses pandas merge. \n    \n    It is similar to using base_shift() we defined earlier from 1 all the way to 100.\n    '''       \n    \n    op = data.merge(\n        data[['Patient', 'Weeks', 'FVC', 'Percent']].rename(\n            {\"Weeks\": \"Base_Week\", \n             \"FVC\": \"Base_FVC\", \n             \"Percent\": \"Base_Percent\"}, axis=1), \n        on='Patient')\n\n    # create week offsets\n    op['Week_Offset'] = op['Weeks'] - op['Base_Week']\n\n    # only take those rows with offset other than 0\n    op = op[op['Week_Offset'] != 0]\n\n    if display:\n        # number of training samples\n        print (\"Number of Samples:{:5} -> {:5}\\nNumber of Columns:{:5} -> {:5}\".format(\n            data.shape[0], op.shape[0], data.shape[1], op.shape[1]))\n    \n    return op.sort_values(by=['Patient', 'Base_Week']).reset_index(drop=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# how would it look?\nop = multi_baseweek_frame(train)\nop.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's write a simple function to preprocess the dataframe to be able to fit to our model:"},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_model_data(data, cat_cols, num_cols, to_drop, cat_method='1h', transform_stats=None,\n                   age_bins=None, train=True, display_stats=True, math=None, factor=False):\n    \n    '''\n    Our pipeline for this notebook. This portion is complex, could be written much more efficiently \n    using simple sklearn tools. I have this bad habit of reinventing the wheel from scratch.\n    \n    Skip this portion. It simply scales, one hot encodes, etc.\n    \n    The main advantage of this function was that it helped me to \n    quickly tweak and see if different preprocessing techniques worked better.\n    '''\n    \n    X = data.copy().reset_index(drop=True)\n    \n    ##########################################################################\n    ##################### NEW FEATURES ADDED COME HERE #######################\n    ##########################################################################\n    \n    \n    if age_bins:    \n        X['binned_age'] = pd.cut(X['Age'], bins=range(0, 101, 100//(age_bins-1))).cat.codes / age_bins\n        to_drop = to_drop + ['Age']   \n        \n    if math: # some simple math based features\n        prod = np.ones(X.shape[0])\n        if 'sin' in math:\n            X['Sin_week'] = X.groupby(\"Patient\")['Weeks'].apply(np.sin)\n            prod = prod * X['Sin_week']\n        if 'cos' in math:\n            X['Cos_week'] = X.groupby(\"Patient\")['Weeks'].apply(np.cos)\n            prod = prod * X['Cos_week']\n        if 'tan' in math:\n            X['Tan_week'] = X.groupby(\"Patient\")['Weeks'].apply(np.tan)\n            prod = prod * X['Cos_week']\n            \n        if len(math) > 1:\n            X['Math_Prod'] = prod\n            \n    if factor:\n        X['factor'] = X['Base_FVC'] / X['Base_Percent']\n        \n    ##########################################################################\n    ##################### ENCODING OF NON NUMERIC DATA #######################\n    ##########################################################################\n    \n    if cat_cols != []:\n    \n        if cat_method == 'ord': # ordinal encoding for tree based models\n            global ordenc\n            if train:\n                from sklearn.preprocessing import OrdinalEncoder\n                ordenc = OrdinalEncoder()\n                X = X.merge(\n                    pd.DataFrame(\n                        ordenc.fit_transform(X[cat_cols]).astype(int),\n                        columns=map(lambda x: x+\"_ord\", cat_cols)),\n                    left_index=True, right_index=True)\n\n            else:\n                X = X.merge(\n                    pd.DataFrame(\n                        ordenc.transform(X[cat_cols]).astype(int),\n                        columns=map(lambda x: x+\"_ord\", cat_cols)),\n                    left_index=True, right_index=True)           \n\n        elif cat_method == '1h': # one hot encoding\n            global onehenc\n            if train:\n                onehenc = OneHotEncoder()\n                X = X.merge(\n                    pd.DataFrame(\n                        onehenc.fit_transform(X[cat_cols]).todense(),\n                        columns=[*np.concatenate(onehenc.categories_)]),\n                    left_index=True, right_index=True)\n\n            else:\n                X = X.merge(\n                    pd.DataFrame(\n                        onehenc.transform(X[cat_cols]).todense(),\n                        columns=[*np.concatenate(onehenc.categories_)]),\n                    left_index=True, right_index=True)\n\n        elif cat_method == 'poly': # polynomial feature encoding\n            global cat_comb\n            if train:\n                cat_comb = np.array(\n                    np.meshgrid(*[X[cat].unique() for cat in cat_cols])\n                ).T.reshape(-1, len(cat_cols))\n\n            for combination in cat_comb:\n                name = \"_\".join(map(str, combination))\n                X[name] = 1\n                for i in range(len(cat_cols)):\n                    X[name] = X[name] & (X[cat_cols[i]] == combination[i]).astype(int)\n\n    # drop the columns after they have been encoded\n    to_drop = to_drop + cat_cols\n                \n    ##########################################################################\n    ######################## SCALING OF NUMERIC DATA #########################\n    ##########################################################################\n    \n    global stats\n    if train:\n        if transform_stats is None:\n            # saving stats for the future\n            stats = X.describe().T\n        else:\n            stats = transform_stats\n\n    # lets scale the numeric columns (We scale it with max possible values)\n    for col in num_cols:\n        \n        if (not train) and (col not in X.columns):\n            continue\n            \n        X[col] = (X[col] - stats.loc[col, 'min']) / (stats.loc[col, 'max'] - stats.loc[col, 'min'])\n        \n    ##########################################################################\n    #################### SPLIT DATAFRAME & DISPLAY STATS #####################\n    ##########################################################################\n\n    global x_cols\n    if train:\n        \n        Y = X['FVC'].dropna()\n        \n        if display_stats:\n        \n            # print out how well our features would do\n            print (X.corr()['FVC'].abs().sort_values(ascending=False)[1:])\n        \n        X = X.drop(to_drop, axis=1)\n        x_cols = X.columns\n        \n        return X, Y\n    \n    else:\n        \n        X = X.drop(to_drop, axis=1, errors='ignore')\n        X = X[x_cols]\n        \n        return X","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's evaluate the above idea to see how it performs:"},{"metadata":{"trusted":true},"cell_type":"code","source":"cat_cols = ['Sex', 'SmokingStatus']\nnum_cols = ['Weeks', 'Week_Offset', 'Base_Week', 'Age', 'Base_FVC', 'Percent', 'Base_Percent']\nto_drop = [\"FVC\", 'Percent', 'Weeks', 'Base_Week', 'Base_Percent']\n\nfolds = 7\nscores = {}\n\n# total patients\ntotal_patients = train.Patient.unique()\nnp.random.shuffle(total_patients)\nval_len = len(total_patients) // folds\n\nX_VAL = base_shift(train, q=0)\nY_VAL = X_VAL['FVC'].dropna()\n\n# most notebooks implement the idea where\n# percent is used in training, & base percent\n# is used for predictions\nX_VAL['Percent'] = X_VAL['Base_Percent']\n\n# creating x_test for predictions simaltanesly\nx_test = sub[['Patient', 'Weeks']].merge(\n    test.rename({\"Weeks\": \"Base_Week\", \n                 \"FVC\": \"Base_FVC\", \n                 \"Percent\": \"Base_Percent\"}, axis=1), \n    on='Patient')\n\n# create week offsets\nx_test['Week_Offset'] = x_test['Weeks'] - x_test['Base_Week']\n# create percent from base_percent\nx_test['Percent'] = x_test['Base_Percent']\n\n# merge data before preprocessing them\nop['Where']     = 'train'\nX_VAL['Where']  = 'val'\nx_test['Where'] = 'test'\ntemp = pd.concat([op, X_VAL, x_test])\n\n# let's create the model data from augmented train\nX, Y = get_model_data(\n    temp, cat_cols=cat_cols, num_cols=num_cols,\n    to_drop=to_drop, display_stats=False, cat_method='1h',\n    math=[],\n)\n\n# get back the data after preprocessing\nx_test = X[X.Where == 'test'].drop('Where', 1).reset_index(drop=True)\nX_VAL = X[X.Where == 'val'].drop(\"Where\", 1).reset_index(drop=True)\nX, Y = X[X.Where == 'train'].drop(\"Where\", 1), Y[X.Where == 'train']\n\npreds = []\n\nfor i in range(folds):\n\n    # do a train/val split for evaluation of the strategy   \n    val_patients = total_patients[(i)*val_len:(i+1)*val_len]\n    train_patients = np.setdiff1d(total_patients, val_patients)\n    \n    x, y = X[X.Patient.isin(train_patients)], Y[X.Patient.isin(train_patients)]\n    x_val, y_val = X_VAL[X_VAL.Patient.isin(val_patients)], Y_VAL[X_VAL.Patient.isin(val_patients)]\n    \n    assert len(np.intersect1d(x_val.Patient.unique(), x.Patient.unique())) == 0\n\n    # how does it perform on train\n    scores['Train'] = scores.get('Train', []) + [cross_val_score(\n        LinearRegression(), x.drop(\"Patient\", 1), y, \n        scoring=l1(70), cv=GroupKFold(5), groups=x.Patient).mean()]\n    \n    # performance on validation data\n    lr = LinearRegression().fit(x.drop(\"Patient\", 1), y)\n    scores['Val'] = scores.get('Val', []) + [-laplace_log_likelihood(\n        y_val, lr.predict(x_val.drop(\"Patient\", 1)), sigma=70\n    ).numpy()]\n    \n    preds.append(lr.predict(x_test.drop(\"Patient\", 1)))\n    \nprint (\"Train Scores Mean: {:.2f} @ {:.2f} Variance\\n\\\nTest Scores Mean : {:.2f} @ {:.2f} Variance\".format(\n        np.mean(scores['Train']), np.std(scores['Train']), \n        np.mean(scores['Val']), np.std(scores['Val'])))\n\npreds = np.stack(preds, 1)\n\n# validation scores\nscores['Val']","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This gives us a higher score than our previous high score of 7.64. Although the train score is lesser, it generalizes better to the test dataset. Let's make predictions with this method and see how it scores on the LB:"},{"metadata":{"trusted":true},"cell_type":"code","source":"# make the predictions\nsub['FVC'] = preds.mean(axis=1)\n\n# altering confidence to get comparable \n# LB score more accurate, lesser confidence\nsub['Confidence'] = 200 + preds.std(axis=1)\n\n# final touches before submission\nfor i in range(len(test)):\n    sub.loc[sub['Patient_Week']==test.Patient[i]+'_'+str(test.Weeks[i]), 'FVC'] = test.FVC[i]\n    sub.loc[sub['Patient_Week']==test.Patient[i]+'_'+str(test.Weeks[i]), 'Confidence'] = 70\n\n# save to csv file\nsub[['Patient_Week', 'FVC', 'Confidence']].to_csv(\"multi_baseweek.csv\", index=False)\n\n# how does it look?\nsub.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Submission Scores:\n- Public LB score: `-6.9577`\n- Private LB score: `-6.8754`\n\nWhich is an improvement of our previous LB high score of `-6.91` @ *375* confidence. \n\n*Although the private LB still remains unbeaten. This probably might be the case since we chose a confidence of 350. Whereas here it's only 200. So the accuracy has greatly improved than before.*\n\nOne of the main reasons why LinearRegression is cool is that it can be used to analyze *how* a set of features combine to produce the result. `Percent` feature is very powerful, let's understand how it works:"},{"metadata":{"trusted":true},"cell_type":"code","source":"temp = base_shift(train, q=0)\n\ntrain_cols = ['Base_Percent', 'Percent', 'Base_FVC']\n\nprint (\"The cross val score is: {}\\n\".format(cross_val_score(\n    LinearRegression(), \n    temp[train_cols], \n    temp['FVC'], \n    scoring=l1(70), cv=GroupKFold(5), \n    groups=train.Patient).mean()))\n\ntemp = LinearRegression().fit(temp[train_cols], temp['FVC'])\nprint (\"The equation is:\\n\", list(zip(train_cols, np.round(temp.coef_, 2))) + [temp.intercept_], sep='')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Can we use the above knowledge for data augmentation?\n\nWe can observe that percent itself does not matter that much. It's relative difference between them that really matters in making FVC predictions. Since Percent is one of the best feature to predict the FVC, by augmenting `Percent` as accurately as possible, we can augment `FVC` as well."},{"metadata":{"trusted":true},"cell_type":"code","source":"def augment_train_cosine(data, n_similar=3, threshold=0.25, display_sample=True):\n    \n    '''\n    - `n_similar` is number of patients we cluster at a time more the cluster, more eratic it gets.\n    - `threshold` is used for as a measure to counter the influence of \"outliers\"\n    '''\n    \n    from sklearn.metrics.pairwise import cosine_similarity\n\n    temp = base_shift(data.copy(), q=0)\n\n    # feature engineering -> simiarity of percentage changes across patients\n    # we use slopes, sex, SmokingStatus and the way the Percent features vary \n    # from Base_Percent as preditive features for clustering patients together.\n    \n    temp['present_minus_past'] = temp.groupby(\"Patient\")['FVC'].transform('diff').fillna(0)\n    temp['Week_diff'] = temp.groupby(\"Patient\")['Weeks'].transform('diff').fillna(0)\n    temp['pms'] = (temp['present_minus_past'] / temp['Week_diff']).replace([np.inf, -np.inf]).fillna(0)\n    temp['pms_min']  = temp.groupby(\"Patient\")['pms'].transform('min')\n    temp['pms_25']   = temp.groupby(\"Patient\")['pms'].transform(lambda x: np.percentile(x, q=25))\n    temp['pms_mean'] = temp.groupby(\"Patient\")['pms'].transform('mean')\n    temp['pms_75']   = temp.groupby(\"Patient\")['pms'].transform(lambda x: np.percentile(x, q=75))\n    temp['pms_max']  = temp.groupby(\"Patient\")['pms'].transform('max')\n    temp['pms_sum']  = temp.groupby(\"Patient\")['pms'].transform('sum')\n    \n    temp['pmb_avg'] = (temp['Percent'] - temp['Base_Percent']).groupby(temp.Patient).transform(\"mean\")\n    temp['p_std']    = temp.groupby(\"Patient\")['Percent'].transform('std')\n\n    temp = temp.merge(\n        pd.concat([\n            temp.groupby(\"Patient\").apply(\n                lambda x: (x['Percent'].values[-1] - x['Percent'].values[0]) / \n                (x['Weeks'].values[-1] - x['Weeks'].values[0])).rename(\"Slope\"),\n            \n            temp.groupby(\"Patient\").apply(\n                lambda x: x['Week_Offset'].iloc[np.argmax(x['pms'])]).rename(\"pmsw_max\"),    \n            temp.groupby(\"Patient\").apply(\n                lambda x: x['Week_Offset'].iloc[np.argmin(x['pms'])]).rename(\"pmsw_min\")\n            \n        ], axis=1, ignore_index=False), on='Patient')\n\n    # we take just the head row for each patient for comparison\n    temp = temp.groupby('Patient').head(1).reset_index(drop=True)\n\n    # another measure of similartiy\n    temp['factor'] = temp['Base_FVC'] / temp['Base_Percent']\n\n    # features to use for similarity clustering\n    train_cols = [\n        # compulsary features\n        'Patient', #'Sex', 'SmokingStatus'\n        \n        # optionally added features for clustering\n        # 'pms_min', 'pms_max', 'p_std',\n        'pms_25', 'pms_mean', 'pms_75', \n        'pms_sum', 'pmsw_min', 'pmsw_max', \n        'pmb_avg', 'Slope', 'factor'\n    ]\n    \n    cat_cols = np.intersect1d(['Sex', 'SmokingStatus'], train_cols)\n\n    temp = temp[train_cols]\n    temp = pd.get_dummies(temp, columns=cat_cols, drop_first=True, prefix='', prefix_sep='')\n\n    # including that patient, find n_similar more patients\n    n_similar += 1\n    groups = (pd.DataFrame(\n        np.argsort(\n            # cosine similarity to get their similarity scores\n            cosine_similarity(temp.drop(\"Patient\", 1), temp.drop('Patient', 1)))\n        [:, -1:-n_similar-1:-1]))\n\n    # cosine similarity is symmetric, so we remove the redundant ones\n    groups = groups[~pd.DataFrame(np.sort(groups.values, axis=1)).duplicated()]\n    \n    # convert the indices to patient ids\n    groups = groups.applymap(lambda x: temp.Patient.to_dict()[x]).apply(list, axis=1).to_dict()\n\n    # the bottle neck of this function\n    aug_data = []\n    for group in tqdm(groups.values(), disable=not display_sample):\n\n        temp = base_shift(train[train.Patient.isin(group)], q=0)\n\n        temp['base_per_diff_from_mean'] = temp['Base_Percent'] - temp['Base_Percent'].unique().mean()\n        temp['Percent_shifted'] = temp['Percent'] - temp['base_per_diff_from_mean']\n\n        temp['base_week_diff_from_mean'] = temp['Base_Week'] - temp['Base_Week'].unique().mean()\n        temp['Week_shifted'] = temp['Weeks'] - temp['base_week_diff_from_mean']\n\n        temp = pd.merge_ordered(\n            temp.drop(['Age', 'Sex', 'SmokingStatus', 'Base_Week', 'Week_Offset'], axis=1),\n\n            # we obtain the mean only for those with samples greater than threshold %\n            # these mean values are fit as such as the expected Percent\n            (temp.groupby(\"Week_shifted\")['Percent_shifted']\n             .agg(['mean', 'count']).query(f'count > {n_similar * threshold}')\n             .drop(\"count\", 1)),\n\n            on='Week_shifted', left_by='Patient'\n        )\n\n        aug_data.append(temp)\n\n    temp = pd.concat(aug_data).reset_index(drop=True)\n\n    # recuring features can simply be padded\n    temp[['base_per_diff_from_mean', 'base_week_diff_from_mean', 'Base_Percent', 'Base_FVC']] = (\n        temp.groupby(\"Patient\")[['base_per_diff_from_mean', 'base_week_diff_from_mean', \n                                 'Base_Percent', 'Base_FVC']].fillna(method='ffill'))\n\n    # get back the weeks from shifted weeks\n    temp['Week_aug'] = temp['Week_shifted'] + temp['base_week_diff_from_mean']\n\n    # For those percent values already present copy them, augment the rest\n    temp['Percent_aug'] = np.where(\n        temp['Percent'].isna(), \n        temp['mean'] + temp['base_per_diff_from_mean'],\n        temp['Percent'])\n\n    # drop/clean those which have neither \n    temp = temp[~temp['Percent_aug'].isna()]\n\n    # fill the FVC values using percent -> FVC corelation\n    test_ids = temp['FVC'].isna()\n    temp.loc[test_ids, 'FVC'] = LinearRegression().fit(\n        temp.loc[~test_ids, ['Base_Percent', 'Percent_aug', 'Base_FVC']], \n        temp.loc[~test_ids, 'FVC']).predict(\n        temp.loc[test_ids, ['Base_Percent', 'Percent_aug', 'Base_FVC']])\n\n    temp = temp.groupby([\"Patient\", 'Week_aug']).mean().reset_index()\n\n    # retain only those columns that we need & rename them to match train\n    temp = temp[['Patient', 'Week_aug', 'Percent_aug', 'FVC']]\n    temp = temp.rename({'Week_aug': 'Weeks', 'Percent_aug': 'Percent'}, axis=1)\n    \n    # some weeks may end up slightly shifted to left or right\n    temp['Weeks'] = temp['Weeks'].astype(int)\n    \n    # add these columns to be able to fit inside multibaseweek pipeline\n    temp = temp.merge(\n        (data[['Patient', 'Sex', 'Age', 'SmokingStatus']]\n         .groupby('Patient').head(1).reset_index(drop=True)), \n        \n        on='Patient')\n    \n    # display how the data has been augmented\n    if display_sample:\n        \n        print (\"Data augmented by factor: {:.2f}x\".format(1 + (\n            temp.shape[0] - data.shape[0]) / data.shape[0]))\n\n        f, ax = plt.subplots(nrows=4, ncols=2, figsize=(20, 20))\n        for i, pat in enumerate(temp.Patient.unique()[:4]):\n\n            ax[i][0].plot(*list(zip(*data[data.Patient == pat][['FVC', 'Weeks']].values))\n                          [::-1], c='g', alpha=0.7)\n\n            ax[i][1].plot(*list(zip(*data[data.Patient == pat][['Percent', 'Weeks']].values))\n                          [::-1], c='g', alpha=0.7)\n\n            ax[i][0].scatter(*list(zip(*temp[temp.Patient == pat][['FVC', 'Weeks']].values))\n                          [::-1], c='r')\n\n            ax[i][1].scatter(*list(zip(*temp[temp.Patient == pat][['Percent', 'Weeks']].values))\n                          [::-1], c='r')\n            \n\n            ax[i][0].set(xlabel='Weeks', ylabel='FVC')\n            ax[i][1].set(xlabel='Weeks', ylabel='Percent')\n            f.suptitle(\"FVC & Percent Augmentation\", y=.9)\n    \n    return temp","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### How does it work?"},{"metadata":{"trusted":true},"cell_type":"code","source":"# default values for n_similar and threshold works fine\n# higher threshold would cause a situation where there are \n# insufficent number of samples\ntemp = augment_train_cosine(train, n_similar=5, threshold=.5)\ntemp.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"*The red dot denotes the augmented data while the green line dentores the original FVC progression.* With increasing thresholds the problem is that data points for a few weeks become large. We can tweak the `train_cols` in the function or add new attributes as a measure of patient similarity.\n\nAnother naive idea of augmenting the data would be interpolate with different methods, let's do that as well and compare the model that performs the best:"},{"metadata":{"trusted":true},"cell_type":"code","source":"def augment_train_naive(\n    data, steps=5, method='index', noise=.25, val_split=0.25, \n    end_pts=[None, None], display_sample=True):\n    \n    '''\n    end_pts -> start and end of augmentation, if None defaults to min/max for that patient\n    '''\n    \n    temp = data[['Patient','Weeks', 'FVC', 'Percent']].merge(\n        \n        ((data.groupby(\"Patient\")['Weeks']\n         .apply(lambda x: pd.Series(\n             np.union1d(np.arange(\n                 end_pts[0] if end_pts[0] else x.min(), \n                 end_pts[1] if end_pts[1] else x.max(), \n                 step=steps), x))\n               ).reset_index(level=0))),\n\n        on=['Patient', 'Weeks'], how='right')\n\n    temp.loc[:, ['FVC', 'Percent']] = (\n        temp.groupby(\"Patient\")[['FVC', 'Percent']]\n        .apply(lambda x: (\n            # interpolate \n            x.interpolate(method=method, limit_direction='both') + \n            \n            # noise factor: Gaussian noice + standard deviation of resp features\n            # (we assume std of percent is scaled up version of std of FVC)\n            (x.std().values * np.random.uniform(-noise, noise, [len(x), 1])))))\n\n    temp = temp.merge(\n        data.groupby(\"Patient\")[['Patient', 'Age', 'Sex', 'SmokingStatus']].head(1),\n        on='Patient')\n    \n    if display_sample:\n        f, ax = plt.subplots(nrows=4, ncols=2, figsize=(20, 20))\n\n        for i, pat in enumerate(temp.Patient.unique()[:4]):\n\n            ax[i][0].plot(*list(zip(*data[data.Patient == pat][['FVC', 'Weeks']].values))\n                [::-1], c='g', alpha=0.7)\n            \n            ax[i][1].plot(*list(zip(*data[data.Patient == pat][['Percent', 'Weeks']].values))\n                          [::-1], c='g', alpha=0.7)\n            \n            ax[i][0].scatter(*list(zip(*temp[temp.Patient == pat][['FVC', 'Weeks']].values))[::-1], c='r')\n            \n            ax[i][1].scatter(*list(zip(*temp[temp.Patient == pat][['Percent', 'Weeks']].values))[::-1], c='r')\n            \n            ax[i][0].set(xlabel='Weeks', ylabel='FVC')\n            ax[i][1].set(xlabel='Weeks', ylabel='Percent')\n            f.suptitle(\"FVC & Percent Augmentation\", y=.9)\n            \n        print (\"Data augmented by factor: {:.2f}x\".format(1 + (\n            temp.shape[0] - data.shape[0]) / data.shape[0]))\n    \n    return temp","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"How does this work?"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Methods to try:\n# 'piecewise_polynomial',, \n# 'pchip', 'akima', 'cubicspline', 'nearest', \n# 'zero', 'slinear', 'quadratic', 'cubic', \n# pad', 'linear'\n\n# for large steps, nearest works best\n# for smaller steps, Linear, Pad works best\n# zero works best for both\n\ntemp = augment_train_naive(\n    data=train, method='akima', steps=15, noise=0\n)\n\ntemp.shape","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's see which aug technique works best. As always we need to create multiple train/val splits to see which one truly works well:"},{"metadata":{"trusted":true},"cell_type":"code","source":"cat_cols = ['Sex', 'SmokingStatus']\nnum_cols = ['Weeks', 'Week_Offset', 'Base_Week', 'Age', 'Base_FVC', 'Percent', 'Base_Percent']\nto_drop = [\"FVC\", 'Percent', 'Weeks', 'Base_Week', 'Age', 'Base_Percent', 'factor']\nmath = []\n\nmulti_method = {}\nmulti_data = {}\nfolds = 7\n\nmethods = ['cubic', 'quadratic', 'cubicspline', \n           'pchip', 'akima', 'nearest', 'zero', \n           'slinear', 'linear', 'SIMILARITY_AUG']\n\n# train patients for CV\ntotal_patients = train.Patient.unique()\nval_len = len(total_patients) // folds\n\nfor method in methods:\n    \n    # create the agumented dataset\n    if method == 'SIMILARITY_AUG':\n        temp = augment_train_cosine(train, n_similar=5, threshold=0.5, display_sample=False)\n    else:\n        temp = augment_train_naive(data=train, method=method, steps=15, noise=0.1, display_sample=False)\n        \n    # create multi base week data\n    temp = multi_baseweek_frame(temp, display=False)\n    # save the augmented dataframe\n    multi_data[method] = temp\n    \n    # creating the data (processed for fitting)\n    X, Y = get_model_data(temp, num_cols=num_cols, cat_cols=cat_cols, to_drop=to_drop, \n        display_stats=False, cat_method='1h', math=math, transform_stats=stats, factor=True)\n        \n    X_VAL = base_shift(train, q=0)\n    Y_VAL = X_VAL['FVC'].dropna()\n    X_VAL['Percent'] = X_VAL['Base_Percent']\n    X_VAL = get_model_data(\n        X_VAL, num_cols=num_cols, cat_cols=cat_cols, factor=True,\n        to_drop=to_drop, train=False, cat_method='1h', math=math)\n    \n    # shuffle the patients for cv\n    np.random.shuffle(total_patients)\n    scores = {}\n    \n    for i in range(folds):\n   \n        val_patients = total_patients[(i)*val_len:(i+1)*val_len]\n        train_patients = np.setdiff1d(total_patients, val_patients)\n\n        x, y = X[X.Patient.isin(train_patients)], Y[X.Patient.isin(train_patients)]\n        x_val, y_val = X_VAL[X_VAL.Patient.isin(val_patients)], Y_VAL[X_VAL.Patient.isin(val_patients)]\n        \n        assert len(np.intersect1d(x_val.Patient.unique(), x.Patient.unique())) == 0\n        \n       # how does it perform on train\n        scores['Train'] = scores.get('Train', []) + [cross_val_score(\n            LinearRegression(), x.drop(\"Patient\", 1), y, \n            scoring=l1(70), cv=GroupKFold(5), groups=x.Patient).mean()]\n\n        # performance on validation data\n        lr = LinearRegression().fit(x.drop(\"Patient\", 1), y)\n        scores['Val'] = scores.get('Val', []) + [-laplace_log_likelihood(\n            y_val, lr.predict(x_val.drop(\"Patient\", 1)), sigma=70\n        ).numpy()]\n        \n    print (\"Method: {}\\nTrain score: {:.2f} @ {:.2f} Variance \\\n    \\nVal Score: {:6.2f} @ {:.2f} Variance\\n{}\\n\".format(\n        method.upper(), np.mean(scores['Train']), np.std(scores['Train']), np.mean(scores['Val']),\n        np.std(scores['Val']), \"=\" * 35\n    ))\n    \n    multi_method[method] = scores","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We need a test score < 7.59 our previous benchmark. \n\nAlmost all our models perform a little better than 7.59 score. We have saved the augmented dataframes. We are going to train models on all these augmented dataframe a aggregate their scores for submission. We can expect the following score on LB @ 70 confidence:"},{"metadata":{"trusted":true},"cell_type":"code","source":"val_scores = [np.mean(multi_method[i]['Val']) for i in methods]\ncut_off = np.percentile(val_scores, 25)\n\nprint(\"At {:.1f} cutoff, the score can be expected to be around: {:.3f}\".format(\n    cut_off, np.mean(list(filter(lambda x: x < cut_off, val_scores)))))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# creating x_test for predictions simaltanesly\nx_test = sub[['Patient', 'Weeks']].merge(\n    test.rename({\"Weeks\": \"Base_Week\", \n                 \"FVC\": \"Base_FVC\", \n                 \"Percent\": \"Base_Percent\"}, axis=1), \n    on='Patient')\n\n# create week offsets\nx_test['Week_Offset'] = x_test['Weeks'] - x_test['Base_Week']\n\n# percent for test\nx_test['Percent'] = x_test['Base_Percent']\n\n# predictions dataframe\npreds = pd.DataFrame(columns=methods)\n\n# train each model on all the saved augmented frames\nfor method in methods:\n    \n    if np.mean(multi_method[method]['Val']) >= cut_off:\n        preds = preds.drop(method, axis=1)\n        continue\n    \n    x, y = get_model_data(\n        multi_data[method], num_cols=num_cols, cat_cols=cat_cols,\n        transform_stats=stats, to_drop=to_drop, train=True, cat_method='1h', \n        display_stats=False, factor=True,\n    )\n\n    lr = LinearRegression().fit(x.drop(\"Patient\", 1), y) \n\n    preds[method] = lr.predict(get_model_data(\n        x_test, cat_cols=cat_cols, num_cols=num_cols, to_drop=to_drop, train=False, factor=True,\n    ).drop(\"Patient\", 1))\n\npreds.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sub['FVC'] = preds.mean(axis=1)\n\n# reduce the confidence to 250\nsub['Confidence'] = 250 + preds.std(1)\n\n# final touches before submission\nfor i in range(len(test)):\n    sub.loc[sub['Patient_Week']==test.Patient[i]+'_'+str(test.Weeks[i]), 'FVC'] = test.FVC[i]\n    sub.loc[sub['Patient_Week']==test.Patient[i]+'_'+str(test.Weeks[i]), 'Confidence'] = 70\n\n# save to csv file\nsub[['Patient_Week', 'FVC', 'Confidence']].to_csv(\"data_aug_submission.csv\", index=False)\n\n# how does it look?\nsub.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Submission Scores:\n- Public LB score: `-6.9089`\n- Private LB score: `-6.8531`\n\nAugmented data is better than score we had obtained by shifting. We say this since we have reduced confidence by 25 and predictions are similar to where we had predicted with 275 confidence.\n\nConfidence values can make or break our model's performance. So let's shift strateges and see some ways how we can implement this feature along with predicting FVC. Some ideas are as follows:\n1. Set confidence values as simply the variance of the the predicted FVC for that person\n2. A secondary model that predicts the confidence scores for a given FVC prediction so as to maximize the metric\n3. Custom training loop with a model predicting FVC and confidence. We use lll loss to train the model.\n\nFor our experimentations let us use new data:"},{"metadata":{"trusted":true},"cell_type":"code","source":"op = multi_baseweek_frame(augment_train_cosine(train, display_sample=False), display=False)\n\n# creating fresh data\nx, y = get_model_data(\n    op, num_cols=num_cols, cat_cols=cat_cols, to_drop=['FVC', 'Percent'], \n    display_stats=False, transform_stats=stats, math=math, factor=True)\n\nx.shape, y.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Implementing Idea 1\nfrom sklearn.model_selection import cross_val_predict\n\ntemp = pd.DataFrame(y)\ntemp['FVC_pred'] = cross_val_predict(\n    LinearRegression(), x.drop([\"Patient\", 'Base_Percent'], 1), \n    y, cv=7, groups=x.Patient)\n\ntemp['Patient'] = x['Patient']\ntemp['Confidence'] = temp.groupby('Patient')['FVC_pred'].transform('std')\ntemp = temp.drop('Patient', axis=1)\n\nprint (laplace_log_likelihood(temp['FVC'], temp['FVC_pred'], sigma=temp['Confidence']).numpy())\ntemp.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Although good, doesn't do all that great. So we are not going to make predictions with this. Let's move onto the next idea, i.e, we predict optimal confidence value for each data in train:"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Implementing Idea 2\ntemp['delta'] = np.minimum((temp['FVC'] - temp['FVC_pred']).abs(), 1000)\n\nfor sigma in range(70, 1000, 10):\n    temp[f'score@{sigma}'] = - (temp['delta'] * np.sqrt(2) / sigma) - np.log(np.sqrt(2) * sigma)\n    \ncols = temp.iloc[:, 4:].columns\n\ntemp['Best_loc'] = temp.iloc[:, 4:].apply(lambda x: x.argmax(), axis=1).astype(int)\ntemp['Best_conf'] = temp['Best_loc'].apply(lambda x: cols[x][6:])\ntemp['Best_score'] = temp.apply(lambda x: x.iloc[x['Best_loc'] + 4], axis=1)\ntemp = temp.iloc[:, [0, 1, 3, 2, -1, -2]].reset_index(drop=True)\n\nprint (\"Laplace Likelihood mean score: {:.3f}\".format(temp['Best_score'].mean()))\ntemp.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"But can we predict this Best_conf ourselves using another model?"},{"metadata":{"trusted":true},"cell_type":"code","source":"# how does it score? RMSE Score\n- cross_val_score(\n    LinearRegression(), \n    \n    # only std as input feature, adding more features doesn't help much\n    temp[['FVC_pred', 'Confidence']], \n\n    temp['Best_conf'].astype(int),\n    scoring='neg_root_mean_squared_error', \n    cv=7\n).mean()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# how much would it benefit our actual score?\n- laplace_log_likelihood(\n    temp['FVC'], temp['FVC_pred'],\n    cross_val_predict(LinearRegression(), temp[['FVC_pred', 'Confidence']], temp['Best_conf'].astype(int))\n).numpy()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# given the data frame with predictions and standard deviation,\n# it outputs the confidence\nconf_model = LinearRegression().fit(\n    temp[['FVC_pred', 'Confidence']],\n    temp[['Best_conf']],\n)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# restoring previous model's predictions\nsub['FVC'] = pd.read_csv(\"multi_baseweek.csv\")['FVC']\n\n# predict and make submission\nsub[['Confidence']] = conf_model.predict(np.stack([\n    sub['FVC'], sub.groupby(\"Patient\")[\"FVC\"].transform(\"std\")], \n        axis=1))\n\n# final touches before submission\nfor i in range(len(test)):\n    sub.loc[sub['Patient_Week']==test.Patient[i]+'_'+str(test.Weeks[i]), 'FVC'] = test.FVC[i]\n    sub.loc[sub['Patient_Week']==test.Patient[i]+'_'+str(test.Weeks[i]), 'Confidence'] = 70\n\nsub[['Patient_Week', 'FVC', 'Confidence']].to_csv(\"data_aug_conf_sub.csv\", index=False)\nsub.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Submission Scores:\n- Public LB score: `-6.9069`\n- Private LB score: `-6.8532`\n\nThis is an improvement from our mutli-weekframe predictions where we used 200 confidence. There's high error in the optimal conf prediction, we could do better if can predict confidence with lesser error.\n\nNow another question arises. Can we do better by simply assigning one conf value per patient?"},{"metadata":{"trusted":true},"cell_type":"code","source":"grouped_scores = pd.DataFrame(columns=['Patient', 'Grp_conf'])\n\nfor i, (pat, group) in enumerate(temp.groupby(x.Patient.reset_index(drop=True))):\n    best_param = [None, np.inf]\n    for sigma in range(70, 1000, 10):\n        curr_score = - laplace_log_likelihood(group['FVC'], group['FVC_pred'], sigma=sigma)\n        if curr_score < best_param[1]:\n            best_param = sigma, curr_score\n    else:\n        grouped_scores.loc[i] = pat, best_param[0]\n    \ntemp = (\n    temp.merge(x.Patient.reset_index(drop=True), left_index=True, right_index=True)\n    .merge(grouped_scores, on=['Patient']).drop('Patient', 1)\n)\n\n# how much would ideally gruped confidence score?\nlaplace_log_likelihood(temp['FVC'], temp['FVC_pred'], temp['Grp_conf'].astype(float)).numpy()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This doesn't perform all that well (optimal CV score is equal to ones we got using a LR in the previous idea), so we drop this one,\n\nLet's now implement the last idea: *Predicting Confidence and FVC in one go.*:\n\nOne way to do that would be to use the [quantile regression](https://scikit-learn.org/stable/auto_examples/ensemble/plot_gradient_boosting_quantile.html) concept where we predict the quuantile intervals and subtract their difference as the confidence score. We are going to use GradientBoostingRegressor for this purpose:\n\n`Note`: GBR by itself didn't prove very helpful. Almost all tree based models performed so poorly. Here are the tricks to make it work:\n1. Choose as max tree depth as possible.\n2. Choose a LR init for the model (Highly imp)\n3. Smaller n_estimators also proved to help CV.\n\nThere are plenty of parameters for `GradientBoostingRegressor`. Therefore to derive on the best set of parameters, we are going to do a gridSearch using the multi_baseweek_frame we had already created, leaving those we know definitely works."},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.ensemble import GradientBoostingRegressor\nfrom sklearn.model_selection import RandomizedSearchCV\n\n# preprocessing pipe essentials\ncat_cols = ['Sex', 'SmokingStatus']\nto_drop = ['FVC', 'Base_Percent', 'Percent']\nnum_cols = ['Weeks', 'Week_Offset', 'Base_Week', 'Age', 'Base_FVC', 'Percent', 'Base_Percent']\ncat_method = 'ord'\nmath = []\nage_bins = 5\n\ngrid_params = {\n    'learning_rate': [5, 1, 0.1, 0.5, 0.075],\n    'subsample': [0.25, 0.5, 0.75, 1.0],\n    'min_samples_split': [2, 0.1, 0.25, 0.5, 0.75, 1.0],\n    'min_samples_leaf': [1, 0.1, 0.25, 0.5, 0.75, 1.0],\n    'min_impurity_decrease': [0.0, 0.1, 0.2],\n    'max_features': ['auto', 'sqrt', 'log2'],\n    'ccp_alpha': [0.0, 0.25, 0.5, 0.75],\n}\n\ngrid = RandomizedSearchCV(\n    \n    GradientBoostingRegressor(\n        random_state=0, loss='lad', \n        n_estimators=50, max_depth=2,\n        criterion='friedman_mse',\n        init=LinearRegression()),\n    \n    param_distributions=grid_params, verbose=1,\n    scoring='neg_root_mean_squared_error', n_iter=300,\n    n_jobs=-1, cv=GroupKFold(n_splits=3), random_state=0)\n\nx, y = get_model_data(\n    op, num_cols=num_cols, cat_cols=cat_cols, age_bins=age_bins, math=math,\n    to_drop=to_drop, display_stats=False, cat_method=cat_method, transform_stats=stats\n)\n\ngrid.fit(x.drop(\"Patient\", 1), y, groups=op.Patient)\nprint (\"Best Score:\", - grid.best_score_)\n\ngrid.best_params_","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Using the above configuration, lets fit and test on the validation data to see how our model would perform."},{"metadata":{"trusted":true},"cell_type":"code","source":"# preprocessing pipe essentials\ncat_cols = ['Sex', 'SmokingStatus']\nto_drop = ['FVC', 'Base_Percent', 'Percent']\nnum_cols = ['Weeks', 'Week_Offset', 'Base_Week', 'Age', 'Base_FVC', 'Percent', 'Base_Percent']\ncat_method = 'ord'\nmath = []\nage_bins = 5\n\n# cross validation \nfolds = 7\ntotal_patients = train.Patient.unique()\nnp.random.shuffle(total_patients)\nval_len = len(total_patients) // folds\n\n# data frame to hold the predictions on train & test data\ntemp = pd.DataFrame()\npreds = pd.DataFrame()\n\n# creating data suitable for model fitting and predictions\nX, Y = get_model_data(\n    op, \n#     base_shift(augment_train_cosine(train, display_sample=False), q=0),\n#     base_shift(train, q=25),\n\n    num_cols=num_cols, cat_cols=cat_cols, age_bins=age_bins, to_drop=to_drop, \n    display_stats=False, cat_method=cat_method, transform_stats=stats, math=math)\n\nX_VAL = base_shift(train, q=0)\nY_VAL = X_VAL['FVC'].dropna()\n\n# percent value as the base percent\nX_VAL['Percent'] = X_VAL['Base_Percent']\n\nX_VAL = get_model_data(\n    X_VAL, num_cols=num_cols, cat_cols=cat_cols, to_drop=to_drop, \n    cat_method=cat_method, train=False, age_bins=age_bins, math=math)\n\n# creating x_test for predictions simaltanesly\nx_test = sub[['Patient', 'Weeks']].merge(\n    test.rename({\"Weeks\": \"Base_Week\", \n                 \"FVC\": \"Base_FVC\", \n                 \"Percent\": \"Base_Percent\"}, axis=1), \n    on='Patient')\n\n# create week offsets\nx_test['Week_Offset'] = x_test['Weeks'] - x_test['Base_Week']\n\n# percent value as the base percent\nx_test['Percent'] = x_test['Base_Percent']\n\nx_test = get_model_data(\n    x_test, cat_cols=cat_cols, num_cols=num_cols, to_drop=to_drop, math=math,\n    train=False, cat_method=cat_method, age_bins=age_bins).drop(\"Patient\", 1)\n\nfor i in range(folds):\n   \n    val_patients = total_patients[(i)*val_len:(i+1)*val_len]\n    train_patients = np.setdiff1d(total_patients, val_patients)\n    \n    assert len(np.intersect1d(val_patients, train_patients)) == 0\n    \n    x, y, = (X[X.Patient.isin(train_patients)].drop(\"Patient\", 1), Y[X.Patient.isin(train_patients)])\n    x_val, y_val = (X_VAL[X_VAL.Patient.isin(val_patients)].drop(\"Patient\", 1), \n                    Y_VAL[X_VAL.Patient.isin(val_patients)])\n        \n    # creating base model, no parameter tweakingparam_distributions\n    model = GradientBoostingRegressor(\n        n_estimators=50, max_depth=2,\n        random_state=0,\n        init=LinearRegression(),\n        criterion='friedman_mse', \n        **grid.best_params_)\n    \n    alpha = 0.75\n\n    model.set_params(loss='quantile', alpha=alpha)\n    model.fit(x.drop([\"Weeks\", 'Base_Week'], 1), y)\n    y_upper = model.predict(x_val.drop([\"Weeks\", 'Base_Week'], 1))\n    y_upper_pred = model.predict(x_test.drop([\"Weeks\", 'Base_Week'], 1))\n\n    model.set_params(loss='quantile', alpha=1-alpha)\n    model.fit(x.drop([\"Weeks\", 'Base_Week'], 1), y)\n    y_lower = model.predict(x_val.drop([\"Weeks\", 'Base_Week'], 1))\n    y_lower_pred = model.predict(x_test.drop([\"Weeks\", 'Base_Week'], 1))\n\n    model.set_params(loss='lad')\n    model.fit(x, y)\n    y_middle = model.predict(x_val)\n    y_middle_pred = model.predict(x_test)\n    \n    print (\"For Fold #{} Val Score: {:.2f} @ 70 Confidence | {:.2f} @ Pred Confidence\".format(\n        i+1, - laplace_log_likelihood(y_middle, y_val), \n        - laplace_log_likelihood(y_middle, y_val, y_upper - y_lower)))\n    \n    temp = temp.append(pd.DataFrame(\n        data=np.stack([y_upper, y_lower, y_middle, y_val], axis=1),\n        columns=['upper', 'lower', 'pred', 'actual']\n    ))\n    \n    preds = preds.append(pd.DataFrame(\n        data=np.stack([y_upper_pred, y_lower_pred, y_middle_pred], axis=1) / folds,\n        columns=['upper', 'lower', 'pred']\n    ))\n    \npreds = preds.groupby(preds.index).sum()\npreds['Confidence'] = preds['upper'] - preds['lower']\n\ntemp['Confidence'] = temp['upper'] - temp['lower']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print (\"\\n|================== Summary ==================|\\n\\\nScore on Total Dataset: {:.3f} @   70 Confidence\\n\\\nScore on Total Dataset: {:.3f} @  225 Confidence\\n\\\nScore on Total Dataset: {:.3f} @ Pred Confidence\".format(\n    -laplace_log_likelihood(temp['actual'], temp['pred'], 70),\n    -laplace_log_likelihood(temp['actual'], temp['pred'], 225),\n    -laplace_log_likelihood(temp['actual'], temp['pred'], temp['Confidence'])\n))\n\nf, ax = plt.subplots(figsize=(40, 40), nrows=4, ncols=2)\nax = ax.ravel()\nfor i, pat in enumerate(np.random.choice(train.Patient.unique(), size=8, replace=False)):\n    (temp.reset_index(drop=True).loc[train.Patient == pat]\n     .drop([\"Confidence\"], 1).plot(ax=ax[i], legend=False))\nf.suptitle(\"Model Predictions\", size=30)\nf.tight_layout(rect=[0, 0.03, 1, 0.95]);","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sub['FVC'] = preds['pred']\nsub['Confidence'] = preds['Confidence']\n\n# final touches before submission\nfor i in range(len(test)):\n    sub.loc[sub['Patient_Week']==test.Patient[i]+'_'+str(test.Weeks[i]), 'FVC'] = test.FVC[i]\n    sub.loc[sub['Patient_Week']==test.Patient[i]+'_'+str(test.Weeks[i]), 'Confidence'] = 70\n\nsub[['Patient_Week', 'FVC', 'Confidence']].to_csv(\"quant_submission.csv\", index=False)\nsub.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"markdown","source":"Submission Scores:\n- Public LB score: `-6.9317`\n- Private LB score: `-6.8600`\n\nMaybe this might work better, if we could really simaltaneously predict confidences instead of doing it seperately. The problem with our linear approaches, in the submission only three weeks are required out of all the 146 weeks of submission. If we are to predict all 146 weeks as close as possible, we might do well with our linear models. Since we are required to predict just 3 of actual submissions (the rest are dropped) we donot have much wiggle space. To get an actual idea, try validating on as little data as possible.\n\n#### Let's combine the two best performing ideas: Data aug + GBR and see how it scores:\n\nThis is the work that you see in my [inference notebook](https://www.kaggle.com/doctorkael/osic-inference) I submitted for this competition. To be able to make predictions on FVC and forecast in one go, we conveniently wrap them up in our custom SciKit Learn class. Methods on how to do this can be found [here](https://scikit-learn.org/stable/developers/develop.html)."},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.base import RegressorMixin, BaseEstimator\nclass GBR(RegressorMixin, BaseEstimator):\n    \n    '''\n    Custom Sckit Learn estimator for making Predictions on FVC and Conf in one go.\n    \n    alpha  -> Quantiles for making predictions\n    umodel -> Makes upper bound preds\n    mmodel -> Makes Middle preds (We optimize this during Grid search)\n    lmodel -> Makes lower bound preds\n    \n    It could be improved by making alpha tunable in grid search. I don't know how \n    this could be done, since cross_val_score accepts only a single pred as op. \n    Maybe we could create a custom scoring function to do this. Let me know if you do ^.^\n    '''\n    \n    def __init__(self, alpha=.75, **params):\n        self.alpha = alpha\n        self.umodel = self._create_model(loss='quantile', q=self.alpha, **params)\n        self.mmodel = self._create_model(loss='lad', **params)\n        self.lmodel = self._create_model(loss='quantile', q=1-self.alpha, **params)\n              \n    def _create_model(self, loss, q=.75, **params):\n        model = GradientBoostingRegressor(\n            init=LinearRegression(),\n            criterion='friedman_mse',\n            n_estimators=50, max_depth=2, \n            loss=loss, alpha=q, **params)\n        \n        return model\n        \n    def fit(self, x, y):\n        '''We fit the same data with all three models.'''\n        \n        self.umodel.fit(x, y)\n        self.mmodel.fit(x, y)\n        self.lmodel.fit(x, y)\n        return self\n    \n    def predict(self, X):\n        '''This function is merely for compatibilty & for providing utility for\n        fine tuning, cross_val_score and as such.'''\n        \n        return self.mmodel.predict(X)\n    \n    def predict_forecast(self, X, return_bounds=False):\n        '''\n        This function make predictions using all three models conveniently wrapped \n        up in one single function. It can be made to produce the bounds or simply the confidence.\n        '''\n        \n        preds = self.mmodel.predict(X)\n        upper = self.umodel.predict(X)\n        lower = self.lmodel.predict(X)\n        \n        if return_bounds:\n            return preds, upper, lower\n        else:\n            return preds, (upper - lower)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Our model is ready. Let's create now cross_val this model on all augmented data to see which aug data performs with least error. The code below may take some time to run:"},{"metadata":{"trusted":true},"cell_type":"code","source":"cat_cols = ['Sex', 'SmokingStatus']\nnum_cols = ['Weeks', 'Week_Offset', 'Base_Week', 'Age', 'Base_FVC', 'Percent', 'Base_Percent']\nto_drop = [\"FVC\", 'Percent', 'Weeks', 'Base_Week', 'Age', 'Base_Percent', 'factor']\nmath = []\n\nmulti_method = {}\nmulti_data = {}\nfolds = 7\n\nmethods = ['cubic', 'quadratic', 'cubicspline', \n           'pchip', 'akima', 'nearest', 'zero', \n           'slinear', 'linear', 'SIMILARITY_AUG']\n\n# train patients for CV\ntotal_patients = train.Patient.unique()\nval_len = len(total_patients) // folds\n\nfor method in methods:\n    \n    # create the agumented dataset\n    if method == 'SIMILARITY_AUG':\n        temp = augment_train_cosine(train, n_similar=3, threshold=0.25, display_sample=False)\n    else:\n        temp = augment_train_naive(data=train, method=method, steps=15, noise=0., display_sample=False)\n        \n    # create multi base week data\n    temp = multi_baseweek_frame(temp, display=False)\n    # save the augmented dataframe\n    multi_data[method] = temp\n    \n    # creating the data (processed for fitting)\n    X, Y = get_model_data(temp, num_cols=num_cols, cat_cols=cat_cols, to_drop=to_drop, \n        display_stats=False, cat_method='1h', math=math, factor=True)\n        \n    X_VAL = base_shift(train, q=0)\n    Y_VAL = X_VAL['FVC'].dropna()\n    X_VAL['Percent'] = X_VAL['Base_Percent']\n    X_VAL = get_model_data(\n        X_VAL, num_cols=num_cols, cat_cols=cat_cols, factor=True,\n        to_drop=to_drop, train=False, cat_method='1h', math=math)\n    \n    # shuffle the patients for cv\n    np.random.shuffle(total_patients)\n    scores = {}\n    \n    for i in range(folds):\n   \n        val_patients = total_patients[(i)*val_len:(i+1)*val_len]\n        train_patients = np.setdiff1d(total_patients, val_patients)\n\n        x, y = X[X.Patient.isin(train_patients)], Y[X.Patient.isin(train_patients)]\n        x_val, y_val = X_VAL[X_VAL.Patient.isin(val_patients)], Y_VAL[X_VAL.Patient.isin(val_patients)]\n        \n        assert len(np.intersect1d(x_val.Patient.unique(), x.Patient.unique())) == 0\n        \n       # how does it perform on train\n        scores['Train'] = scores.get('Train', []) + [cross_val_score(\n            GBR(), x.drop(\"Patient\", 1), y, \n            scoring=l1(70), cv=GroupKFold(5), groups=x.Patient).mean()]\n\n        # fit to measure model's performance\n        lr = GBR().fit(x.drop(\"Patient\", 1), y) \n        temp = lr.predict_forecast(x_val.drop(\"Patient\", 1))\n        temp = pd.DataFrame(np.stack(temp, 1), columns=['pred', 'conf'])\n        temp['actual'] = y_val.reset_index(drop=True)\n        \n        # performance on validation data\n        scores['Val'] = scores.get('Val', []) + [-laplace_log_likelihood(\n            temp['actual'], temp['pred'], 70\n        ).numpy()]\n        \n        # performance with confidence        \n        scores['ValC'] = scores.get('ValC', []) + [-laplace_log_likelihood(\n            temp['actual'], temp['pred'], temp['conf']\n        ).numpy()]\n        \n        # Worst Performing mean scores\n        scores['ValW'] = scores.get('ValW', []) + [temp.apply(lambda x: -laplace_log_likelihood(\n            x['actual'], x['pred'], x['conf']).numpy(), 1).nlargest(25).mean()]\n        \n    print (\"Method: {}\\nTrain score: {:5.2f} @ {:.2f} Variance \\\n    \\nVal Score: {:7.2f} @ {:.2f} Variance\\\n    \\nValC Score: {:6.2f} @ {:.2f} Variance\\\n    \\nWorst Score: {:5.2f} @ {:.2f} Variance\\n{}\\n\".format(\n        method.upper(), np.mean(scores['Train']), np.std(scores['Train']), \n        np.mean(scores['Val']), np.std(scores['Val']), np.mean(scores['ValC']), \n        np.std(scores['ValC']), np.mean(scores['ValW']), np.std(scores['ValW']), \"=\" * 35\n    ))\n    \n    multi_method[method] = scores","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We aim to find those augmented data which yeilds the lowest possible error in validation data. Previously we aimed to increase the score, now we aim to decrease the worst. We also increase the threshold than before and retain 50% of top performing data.\n\n*The rationale for this decision is that in the end, out of all our predictions, only three week FVC are computed for the score.*"},{"metadata":{"trusted":true},"cell_type":"code","source":"val_scores = [np.mean(multi_method[i]['ValW']) for i in methods]\ncut_off = np.percentile(val_scores, 50)\n        \nprint(\"At {:.1f} cutoff, the Worst Score would be {:.3f} @ Pred Conf\".format(\n    cut_off, np.mean(list(filter(lambda x: x < cut_off, val_scores)))))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# creating x_test for predictions simaltanesly\nx_test = sub[['Patient', 'Weeks']].merge(\n    test.rename({\"Weeks\": \"Base_Week\", \n                 \"FVC\": \"Base_FVC\", \n                 \"Percent\": \"Base_Percent\"}, axis=1), \n    on='Patient')\n\n# create week offsets\nx_test['Week_Offset'] = x_test['Weeks'] - x_test['Base_Week']\n\n# percent for test\nx_test['Percent'] = x_test['Base_Percent']\n\n# predictions dataframe\npreds = {i: None for i in methods}\n\n# train each model on all the saved augmented frames\nfor method in methods:\n    \n    if np.mean(multi_method[method]['ValW']) >= cut_off:\n        preds.pop(method)\n        continue\n    \n    x, y = get_model_data(\n        multi_data[method], num_cols=num_cols, cat_cols=cat_cols, \n        to_drop=to_drop, train=True, cat_method='1h', \n        display_stats=False, factor=True,\n    )\n\n    lr = GBR().fit(x.drop(\"Patient\", 1), y) \n\n    preds[method] = lr.predict_forecast(get_model_data(\n        x_test, cat_cols=cat_cols, num_cols=num_cols, to_drop=to_drop, train=False, factor=True,\n    ).drop(\"Patient\", 1), return_bounds=True)\n\npreds.keys()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"temp = {}\nfor i in multi_method:\n    if i in preds.keys():\n        temp['Val'] = temp.get('Val', []) + [np.mean(multi_method[i]['Val'])]\n        temp['ValC'] = temp.get('ValC', []) + [np.mean(multi_method[i]['ValC'])]\n        \nprint (\"At same cutoff:\\n\\nBest val score: {:6.3f} @   70 Conf\\n\\\nBest Val score: {:6.3f} @ Pred Conf\".format(np.mean(temp['Val']), np.mean(temp['ValC'])))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"##### A note on what is done in snippet below: \n\nWe produce the bounds (lower and upper) using models trained on all the augmented data. For each week, we take the minimum (for lower bound) and maximum (for upper bound) to arrive at the actual confidence in that particular predictions."},{"metadata":{"trusted":true},"cell_type":"code","source":"mean = np.mean(np.stack(pd.DataFrame(preds).iloc[0].values), 0)\nupper = np.max(np.stack(pd.DataFrame(preds).iloc[1].values), 0)\nlower = np.min(np.stack(pd.DataFrame(preds).iloc[2].values), 0)\n\nsub['FVC'] = mean\nsub['Confidence'] = upper - lower\n\n# save to csv file\nsub[['Patient_Week', 'FVC', 'Confidence']].to_csv(\"submission.csv\", index=False)\n\nsub.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"markdown","source":"Submission Scores:\n- Public LB score: `-6.8897`\n- Private LB score: `-6.8479`\n\n*This was the submission file that I submitted*\n\n### Some more experiments:\n\n*My experiments after this point weren't implemented or tested by making submissions. It might have performed better or could have performed worse. But I document them here neverthless.*\n\n*Some were used in creating features such as `factor`. Some were used in creating the cosine augmented data. Most were unimplemented ideas. :(*\n\nThe next idea is to use the images of scans to learn the slope before hand & try to predict the FVC. Let's get started:"},{"metadata":{"trusted":true},"cell_type":"code","source":"temp = train.copy() \ntemp['present_minus_past'] = temp.groupby(\"Patient\")['FVC'].transform('diff').fillna(0)\ntemp['Week_diff'] = temp.groupby(\"Patient\")['Weeks'].transform('diff').fillna(0)\ntemp['pms'] = (temp['present_minus_past'] / temp['Week_diff']).replace([np.inf, -np.inf]).fillna(0)\n\n# features we are trying to predict from the image\ntemp['min_pms'] = temp.groupby(\"Patient\")['pms'].transform('min')\ntemp['avg_pms'] = temp.groupby(\"Patient\")['pms'].transform('mean')\ntemp['max_pms'] = temp.groupby(\"Patient\")['pms'].transform('max') \n\n# # these features can be engineered\ntemp['cum_pms'] = (temp['avg_pms'] * temp['Week_diff']).groupby(temp.Patient).transform('cumsum')\ntemp['cum_min_pms'] = (temp['min_pms'] * temp['Week_diff']).groupby(temp.Patient).transform('cumsum')\ntemp['cum_max_pms'] = (temp['max_pms'] * temp['Week_diff']).groupby(temp.Patient).transform('cumsum')\ntemp['Base_FVC'] = temp.groupby(\"Patient\")['FVC'].transform(\"first\")\ntemp['Base_Percent'] = temp.groupby(\"Patient\")['Percent'].transform(\"first\")\n\n# how does it look?\n(temp[temp.Patient == np.random.choice(temp.Patient)]\n [['Base_FVC', 'Weeks', 'min_pms', 'avg_pms', 'max_pms', \n   'cum_min_pms', 'cum_pms', 'cum_max_pms', 'FVC']])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# if successful, we would get a very good score\nprint (\"The cross val score is: {:.2f}\".format(\n    cross_val_score(\n        LinearRegression(),\n        temp[['Base_FVC', 'cum_min_pms', 'cum_max_pms', 'cum_pms']], \n        temp['FVC'], scoring=l1(70), groups=temp.Patient, \n        cv=GroupKFold(7)).mean()))\n\nlr = LinearRegression().fit(temp[['Base_FVC', 'cum_pms']], temp['FVC'])\nprint (\"The equation is:\", list(zip(['Base_FVC', 'cum_pms'], np.round(lr.coef_, 2))) + [lr.intercept_], sep='\\n')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Remember the previous output we had? This is eerily similar to that one. **Guess we have figured out what Percent value is all about.** Although the performance is poor. The previous equation we had is:\n\n```\nThe cross val score is: 5.260435581207275\n\nThe equation is:\n[('Base_Percent', -31.54), ('Percent', 32.36), ('Base_FVC', 0.97), 14.667705269276212]\n```\n\n$ => (Percent * 32.36) - (Base\\_Percent * 31.54) = cum\\_pms * 0.28 $ \n\n$ =>(Percent\\_cum\\_sum * 32.36) + (Base\\_Percent * 0.82) = cum\\_pms * 0.28 $ \n\nLet's verify this equation by predicting percent ourselves:"},{"metadata":{"trusted":true},"cell_type":"code","source":"temp['Percent_pred'] = ((temp['cum_pms'] * 0.28) + (temp['Base_Percent'] * 31.54)) / 32.36\n\n# root mean squared error between them\nnp.sqrt(np.mean((temp['Percent'] - temp['Percent_pred']) ** 2))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This is tremendously accurate, although lets see if we can do better. Some more analysis and probing:"},{"metadata":{"trusted":true},"cell_type":"code","source":"temp = base_shift(train, q=0)\n# preset minus past\ntemp['pms_percent'] = temp.groupby(\"Patient\")['Percent'].transform('diff').fillna(0)\ntemp['pms_FVC'] = temp.groupby('Patient')['FVC'].transform(\"diff\").fillna(0)\ntemp['pms_ratio'] = (temp['pms_FVC'] / temp['pms_percent']).fillna(0)\n\n# present minus base\ntemp['pmb_percent'] = temp.groupby(\"Patient\").apply(\n    lambda x: x['Percent'] - x['Base_Percent']).reset_index(level=0)[0]\n\ntemp['pmb_FVC'] = temp.groupby(\"Patient\").apply(\n    lambda x: x['FVC'] - x['Base_FVC']).reset_index(level=0)[0]\n\ntemp['pmb_ratio'] = (temp['pmb_FVC'] / temp['pmb_percent']).fillna(0)\n\ntemp['fvc_percent_ratio'] = temp['FVC'] / temp['Percent']\n\ntemp[['pms_ratio', 'pmb_ratio', 'fvc_percent_ratio']].sample(5)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We observe that the all three ratios are identical for all three cases! Which means one thing:"},{"metadata":{"trusted":true},"cell_type":"code","source":"temp = base_shift(train, q=0)\n\ntemp['factor'] = temp['Base_FVC'] / temp['Base_Percent']\n\n# the FVC Equation from Percent\ntemp['FVC_pred'] = temp['Percent'] * temp['factor']\n\n# Total error in predictions\n(temp['FVC_pred'] - temp['FVC']).sum()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The FVC can be derived from Percent through the following formula:\n\n$ FVC => Percent * \\dfrac{Base\\_FVC}{Base\\_Percent} $\n\n*It's only at this point that I discovered that Percent is simply a **factor** away from the FVC itself.*\n\nLet's check for some corelation for the `factor` feature. We observe that factor by itself has a good corelation to FVC."},{"metadata":{"trusted":true},"cell_type":"code","source":"(\n    temp.drop(['Base_Week', 'Base_FVC', 'FVC_pred', 'Base_Percent', 'Week_Offset'], 1)\n    .corr()\n)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"How well would a LR with polynomial features preprocessing step do?"},{"metadata":{"trusted":true},"cell_type":"code","source":"# prepare data with LR model\ntemp = temp.drop(['FVC_pred'], 1)\n\nx = pd.get_dummies(\n    temp, columns=['Sex', 'SmokingStatus'], drop_first=True,\n    prefix='', prefix_sep='')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.preprocessing import PolynomialFeatures\nfrom sklearn.decomposition import PCA\n\nmodel = Pipeline([ \n    ('poly', PolynomialFeatures(degree=2)), \n    (\"pca\", PCA(n_components=10)),\n    (\"lr\", LinearRegression())])\n\nprint (\"This Pipeline would produce {} features & scores: {:.2f}, \\\nwhereas a vanilla LR scores: {:.2f}\".format(\n    model.steps[1][1].fit_transform(\n        model.steps[0][1].fit_transform(\n        x.drop([\"Patient\", \"FVC\", 'Percent', 'Base_FVC', 'Weeks', 'Base_Week'], 1))).shape[1],\n    \n    cross_val_score(\n    model, x.drop([\"Patient\", 'FVC', 'Percent', 'Base_FVC', 'Weeks', 'Base_Week'], 1), \n    x['FVC'], cv=GroupKFold(5), groups=x.Patient,\n    scoring=l1(70)).mean(),\n    \n    cross_val_score(\n    LinearRegression(), x.drop([\"Patient\", 'FVC', 'Percent', 'Base_Week', 'Weeks'], 1), \n    x['FVC'], cv=GroupKFold(5), groups=x.Patient,\n    scoring=l1(70)).mean()  \n))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Not much of an improvement. So we drop this idea. Let's resume our previous idea to create features such as `pms_avg`, `pmb_avg`, etc:"},{"metadata":{"trusted":true},"cell_type":"code","source":"from scipy.stats import linregress\ntemp['pmb_avg'] = (temp['Percent'] - temp['Base_Percent']).groupby(temp.Patient).transform(\"mean\")\n\ntemp['pms'] = (\n    temp.groupby(\"Patient\")['Percent'].diff().fillna(0) / \n    temp.groupby(\"Patient\")['Weeks'].diff().fillna(0)\n).replace([-np.inf, np.inf], np.nan).fillna(0)\n\ntemp['pms_avg'] = temp.groupby(\"Patient\")['pms'].transform('mean')\n\ntemp = temp.groupby(\"Patient\").head(1).reset_index(drop=True)\n\ntemp = temp.merge(pd.concat([\n    \n    train.groupby(\"Patient\").apply(\n        lambda x: pd.Series(linregress(x['Weeks'], x['FVC'])[:2], \n                            index=['slope', 'I_FVC'])),\n    \n    train.groupby(\"Patient\")['FVC'].mean().rename(\"FVC_mean\"),\n    train.groupby(\"Patient\")['FVC'].min().rename(\"FVC_min\"),\n    train.groupby(\"Patient\")['FVC'].max().rename(\"FVC_max\"),\n    \n    ], axis=1),\n    on='Patient'\n)\n\n# drop unneccessary features\ntemp = temp.drop(['pms', 'FVC', 'Weeks', 'Percent', 'Week_Offset'], 1)\n\ntemp.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"A bit of Eda on this new feature: Factor:"},{"metadata":{"trusted":true},"cell_type":"code","source":"# unique patiets, unique factors\ntemp.Patient.nunique(), temp.factor.nunique()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"temp.factor.plot(kind='kde', color='yellow', secondary_y=True, title='Factor Distribution Plot')\ntemp.factor.plot(kind='hist', figsize=(15, 5), bins=15, color='r', grid=True);","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"ages = pd.cut(temp.Age.clip(0, 100), bins=[0, 60, 65, 70, 75, 80, 85, 90, 100]).cat.codes\nplt.figure(figsize=(15, 5))\nsc = plt.scatter(temp.index, temp.factor, c=temp.Sex.map({\"Male\": 1, \"Female\": 0}), s=ages*15)\nplt.legend(handles=sc.legend_elements()[0], labels=['Female', 'Male'])\ntemp.groupby([\"Sex\"])['factor'].agg(['min', 'mean', 'max', 'std', 'count'])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"smkunique = temp.SmokingStatus.unique()[::-1]\nplt.figure(figsize=(15, 5))\nsc = plt.scatter(temp.index, temp.factor, alpha=.8, s=ages*15,\n                 c=temp.SmokingStatus.map(dict(zip(smkunique, range(3)))))\n\nplt.legend(handles=sc.legend_elements()[0], labels=smkunique.tolist())\ntemp.groupby([\"SmokingStatus\"])['factor'].agg(['min', 'mean', 'max', 'std', 'count'])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"It is pretty apparant from the above graph, that factor is linearly seperable with Sex, plus Younger people are given higher percent values compared to older people. \n\nThere's no observable relationship among the different SmokingGroups."},{"metadata":{"trusted":true},"cell_type":"code","source":"# saving values for future use\nbase = temp.copy()\nbase['Recovering'] = tf.nn.sigmoid(base['slope']).numpy()\n\ncorr = base.drop(['FVC_min', 'FVC_max'], 1).corr()\n\n# this is really cool way to viz corelation\n# its even better than plt.imshow \ncorr.style.background_gradient(cmap='coolwarm')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"From what we observe, factor can be sorted in ascending to descending order as follows:\n\n1. Female + Older -> Lowest factor\n2. Female + Younger -> Moderate low factor\n3. Male + Older -> Moderate high factor\n4. Male + Younger -> Highest factor\n\nFrom the above correlation map, we also observe that another feature is involved in this equation: `Base_FVC`."},{"metadata":{"trusted":true},"cell_type":"code","source":"ages = pd.cut(temp.Age.clip(0, 100), bins=[0, 60, 65, 70, 75, 100]).cat.codes\nplt.figure(figsize=(15, 5))\nsc = plt.scatter(temp.Base_FVC, temp.factor, c=temp.Sex.map({\"Male\": 1, \"Female\": 0}), s=ages*15)\nplt.legend(handles=sc.legend_elements()[0], labels=['Female', 'Male'])\nplt.title(\"Complete factors involved behind feature: `Factor`\");","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The graph says it all. Greater the Base_FVC, greater the factor. Age, Sex are furthermore used for computing factor."},{"metadata":{"trusted":true},"cell_type":"code","source":"temp.sample(10)[['Base_Percent', 'Base_FVC', 'Sex', 'Age', 'factor', 'slope']].sort_values(\n    ['Sex', 'Age'], ascending=[True, False])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# use pmb_avg & pms_avg when using base_shift\ntemp = pd.get_dummies(base_shift(train, q=0).merge(\n        base.iloc[:, [0, 7, 8, 9, 10, 11, 12, 13, 14]],\n        on='Patient'), columns=['Sex', 'SmokingStatus'], \n               drop_first=True, prefix='', prefix_sep='')\n\n# we do a cumsum here since we can then subtract the cumsum\n# from base FVC when needed\ntemp['pms_avg_cum'] = (\n    temp.groupby(\"Patient\")['Weeks'].transform(\"diff\").fillna(0) * \n    temp['pms_avg']\n).groupby(temp.Patient).cumsum()\n\ntemp['Recovering'] = temp['slope'] > 0\n\ntemp.columns","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's see if we can predict these features using a LR model and then try to make predictions using it (Essentially chaining predictions):"},{"metadata":{"trusted":true},"cell_type":"code","source":"to_drop = [\"Patient\", 'Percent', 'FVC', 'Weeks', 'Base_Week',\n           'pmb_avg', 'pms_avg_cum', 'slope', 'I_FVC', 'FVC_mean', 'FVC_min', \n           'FVC_max', 'pms_avg']\n\nprint (\"Baseline score to beat: {:.2f}\\n\".format(\n     cross_val_score(\n        LinearRegression(),\n        temp.drop(to_drop + ['Male', 'Ex-smoker', 'Never smoked', 'Age'], 1),\n        temp['FVC'],\n        groups=temp.Patient,\n        scoring=l1(70),\n        cv=GroupKFold(5)).mean()\n))\n\nfor no_drop in [['pmb_avg'], ['pms_avg'], ['pms_avg_cum'], \n                ['slope', 'I_FVC'], ['FVC_mean'], ['FVC_min', 'FVC_max']]:\n    \n    if no_drop in [['FVC_ratio'], ['pms_avg_cum'], ['slope'], ['pms_avg']]:\n        To_drop = to_drop + ['Male', 'Ex-smoker', 'Never smoked', 'Age']\n    else:\n        To_drop = to_drop\n\n    print (\"The score with {:20} as features is: {:.2f} (Actual) | {:.2f} (Pred) | {:.2f} (Avg)\"\n           .format(\n               ', '.join(no_drop), \n        \n            cross_val_score(\n            LinearRegression(),\n            temp.drop(np.setdiff1d(To_drop, no_drop), 1),\n            temp['FVC'],\n            groups=temp.Patient,\n            scoring=l1(70),\n            cv=GroupKFold(5)).mean(),\n\n            cross_val_score(\n            LinearRegression(),\n            np.hstack([temp.drop(To_drop, 1), cross_val_predict(\n                LinearRegression(),\n                temp[['Base_FVC', 'Base_Percent', \n                      'Week_Offset', 'factor', \n                      'Age', \"Male\"]],\n                temp[no_drop],\n                cv=GroupKFold(5),\n                groups=temp.Patient,)]),\n            temp['FVC'],\n            groups=train.Patient,\n            scoring=l1(70),\n            cv=GroupKFold(5)).mean(),\n               \n            cross_val_score(\n            LinearRegression(),\n            pd.concat([\n                temp.drop(To_drop, 1), \n                temp.groupby(['Male'])[no_drop].transform('mean')\n            ], axis=1),\n            temp['FVC'],\n            groups=train.Patient,\n            scoring=l1(70),\n            cv=GroupKFold(5)).mean()\n    ))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"markdown","source":"Although they are successful in beating our baseline score, they are still unsuccessful in beating our augment-multibaseweek pipeline. Let's try to predict these features from the image scans to the best we can."},{"metadata":{"trusted":true},"cell_type":"code","source":"base.sample(5).drop(\"Patient\", 1)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Of all these features, the most important feature that makes sense to predict from the scans is `pmb_avg` & `slope`. Let's do a bit of analysis on the data we hadn't analysed before we proceed:"},{"metadata":{"trusted":true},"cell_type":"code","source":"base['pmb_avg'].plot(kind='hist', figsize=(15, 5), xticks=range(-150, 150, 15), bins=25, color='r')\nbase['pmb_avg'].plot(kind='kde', secondary_y=True, title='Avg_pmb Plot');","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"f, ax = plt.subplots(nrows=5, ncols=2, figsize=(20, 25), sharex=False, sharey=False)\n\nfeatures = ['FVC_mean', 'FVC_min', 'FVC_max', 'slope', 'pmb_avg']\nfor i, feat in enumerate(features):\n    base[feat][base.Sex == 'Male'].plot(kind='hist', bins=50, ax=ax[i][0], color='red')\n    base[feat][base.Sex == 'Male'].plot(kind='kde', secondary_y=True, ax=ax[i][0])\n    ax[i][0].set(title=feat+'_male')\n    \n    base[feat][base.Sex == 'Female'].plot(kind='hist', bins=50, ax=ax[i][1], color='blue')\n    base[feat][base.Sex == 'Female'].plot(kind='kde', secondary_y=True, ax=ax[i][1])\n    ax[i][1].set(title=feat+'_female')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can observe that *Females* distribution are much more scattered than it is for Males due to lesser samples available. However the distribution themselves are the same."},{"metadata":{"trusted":true},"cell_type":"code","source":"f, ax = plt.subplots(nrows=5, ncols=3, figsize=(22, 25), sharex=False, sharey=False)\n\nfor i, feat in enumerate(features):\n    base[feat][base.SmokingStatus == 'Currently smokes'].plot(kind='kde', secondary_y=True, ax=ax[i][0])\n    base[feat][base.SmokingStatus == 'Currently smokes'].plot(kind='hist', bins=50, ax=ax[i][0], color='red')\n    ax[i][0].set(title=feat+'_Currently_smoking')\n    \n    base[feat][base.SmokingStatus == 'Never smoked'].plot(kind='kde', secondary_y=True, ax=ax[i][2])\n    base[feat][base.SmokingStatus == 'Never smoked'].plot(kind='hist', bins=50, ax=ax[i][2], color='green')\n    ax[i][2].set(title=feat+'_Never_smoked')\n    \n    base[feat][base.SmokingStatus == 'Ex-smoker'].plot(kind='kde', secondary_y=True, ax=ax[i][1])\n    base[feat][base.SmokingStatus == 'Ex-smoker'].plot(kind='hist', bins=50, ax=ax[i][1], color='yellow')\n    ax[i][1].set(title=feat+'_Ex-smoker')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Let's start with the Scans now:\n\nLets start with the images we have been provided with now. How many slices per patient exists?"},{"metadata":{"trusted":true},"cell_type":"code","source":"temp = {}\nfor pat, _ in map(lambda x: x.split(\"/\")[-2:], train_files):\n    temp[pat] = temp.get(pat, 0) + 1\n    \nplt.figure(figsize=(20, 5))\nplt.hist(temp.values(), bins=50)\nplt.xticks(range(0, 1051, 50));","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"markdown","source":"ID's without an ImagePositionPatient attribute: ID00128637202219474716089, ID00132637202222178761324, ID00026637202179561894768\n\nID's that require GDCM: ID00011637202177653955184, ID00052637202186188008618"},{"metadata":{"trusted":true},"cell_type":"code","source":"BAD_IDS = ['ID00011637202177653955184', 'ID00052637202186188008618']","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Let's write a function to extract the scan images in order. \n\nWe want to be able to stride and load them when required. Say just every 5th image, every 10th image and so on. So we use a parameter called `max_size` to load specified strides of images. Although some or maybe a lot of data might lost in this process, it helps to have all scans of a particular size. We chose size = 10 as default since the minimum scans for a patient is 12."},{"metadata":{"trusted":true},"cell_type":"code","source":"def extract_patient(patient_loc, max_size=10):    \n    import glob\n    \n    images, positions, instances = [], [], []\n    for dcm_file in glob.glob(patient_loc+\"/*\"):\n        temp = pydicom.dcmread(dcm_file)\n        \n        if 'ImagePositionPatient' in temp:\n                positions.append(temp.ImagePositionPatient[2])\n                \n        instances.append(temp.InstanceNumber)\n        images.append(temp.pixel_array)\n            \n    # for those images where position exists for all images\n    if len(positions) == len(images):\n        patients = np.argsort(positions)\n    \n    # certain images don't have this attrib, and we have to use\n    # the Instance number for ordering the scan images\n    else:\n        patients = np.argsort(instances)\n        \n    if max_size == None: # max_size = None -> extract all images\n        max_size = len(images)\n        \n    # here we stride and load the images, dropping those in between\n    patients = patients[np.linspace(0, len(images)-1, max_size).astype(int)]\n    return  np.stack(images)[patients]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"images = extract_patient(\n    f'{main_dir}/train/ID00052637202186188008618', # a bad Image, check if GDCM works\n    max_size=10)\n\nimages.shape, images.min(), images.max()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"f, ax = plt.subplots(nrows=2, ncols=5, figsize=(15, 5))\nax = ax.ravel()\nfor i, img in enumerate(images):\n    ax[i].imshow(img, cmap='bone')\n    ax[i].axis('off')\n    ax[i].set(title=str(i + 1)+'th Image')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's visualize their std, min, mean and max. Maybe this could simply be fed to our model in making its predictions."},{"metadata":{"trusted":true},"cell_type":"code","source":"f, ax = plt.subplots(ncols=4, figsize=(20, 5))\nfor i, (name, img) in enumerate([\n    (\"Min\", images.min(0)), (\"Mean\", images.mean(0)), \n    (\"Max\", images.max(0)), (\"Std\", images.std(0))]):\n\n    ax[i].imshow(img, cmap='bone')\n    ax[i].axis('off')\n    ax[i].set(title=name)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's now create a function to visualize the scans as an animation. Code copied in part from [here](https://www.kaggle.com/gunesevitan/osic-pulmonary-fibrosis-progression-eda):"},{"metadata":{"trusted":true},"cell_type":"code","source":"def display_scan(images):\n    '''\n    We alter and modify gunesevitan's function to be able to display the animation\n    according to the number of slices we have in the scan.\n    '''\n    import matplotlib.animation as animation\n\n    fig, ax = plt.subplots(figsize=(7, 7))\n\n    ims = []\n    for i in images:\n        im = ax.imshow(i, animated=True, cmap=plt.cm.bone)\n        ax.axis('off')\n        ims.append([im])\n\n    ani = animation.ArtistAnimation(fig, ims, interval=1500/len(images), repeat_delay=1000)\n    return ani.to_html5_video()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%capture\nani = display_scan(images)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from IPython.display import HTML\nHTML(ani)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The code below can be used for segmenting the scans. \n\n*I am total newb so I may be totally incorrect in doing this. But doing this helped in improving my embeddings loss a bit more than feeding in the scans after a simple scaling. Maybe my segmentation can be still improved by following all the procdures from [here](https://www.raddq.com/dicom-processing-segmentation-visualization-in-python/). I didn't have much time to play around.*"},{"metadata":{"trusted":true},"cell_type":"code","source":"def segment_scans(images, n_clusters=2):\n    '''Function to segment the Lungs from background'''\n    from sklearn.cluster import KMeans\n    \n    # we do this so that, KMeans only sees the lung area\n    q, (timesteps, shape) = .25, images.shape[:2]\n    first, last = int(shape * q), int(shape * (1 - q))\n    \n    # find out the threshold for that scan using KMeans\n    km = KMeans(n_clusters=n_clusters)\n    cc = km.fit(images[:, first:last, first:last].reshape(timesteps, -1)).cluster_centers_\n    threshold = np.mean(cc)\n    \n    return np.where(images < threshold, 1., 0.)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%capture\nsegmented_images = segment_scans(images / 2048)\nani = display_scan(segmented_images)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"How the same image looks after segmentation:"},{"metadata":{"trusted":true},"cell_type":"code","source":"print (segmented_images.min(), segmented_images.max(), np.unique(segmented_images))\nHTML(ani)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"markdown","source":"Although unsegmented images intially perform good, our segmented images model is able to better capture the scans as embeddings. *Perhaps we need a more powerful model with unsegmented images.*\n\nFrom these images as we had mentioned earlier, we try to predict two values:\n1. slope of FVC of the Patient\n2. Recovering or not\n\nFor this purpose we create an **Convolutional Autoencoder.**"},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_model_small(timesteps=10, imgsize=256, segment=True):\n    \n    seq = tf.keras.models.Sequential()\n    seq.add(tf.keras.layers.Input(shape=(timesteps, imgsize, imgsize, 1)))\n    seq.add(tf.keras.layers.TimeDistributed(\n        tf.keras.layers.Conv2D(64, (11, 11), strides=2, padding=\"same\")))\n    seq.add(tf.keras.layers.LayerNormalization())\n    seq.add(tf.keras.layers.TimeDistributed(\n        tf.keras.layers.Conv2D(32, (7, 7), strides=2, padding=\"same\")))\n    seq.add(tf.keras.layers.LayerNormalization())\n    seq.add(tf.keras.layers.TimeDistributed(\n        tf.keras.layers.Conv2D(16, (5, 5), strides=2, padding=\"same\")))\n    seq.add(tf.keras.layers.LayerNormalization())\n    \n    # # # # #\n    seq.add(tf.keras.layers.ConvLSTM2D(16, (3, 3), padding=\"same\", return_sequences=True))\n    seq.add(tf.keras.layers.LayerNormalization())\n    seq.add(tf.keras.layers.ConvLSTM2D(8, (3, 3), padding=\"same\", return_sequences=True, name='enc_op'))\n    seq.add(tf.keras.layers.LayerNormalization())\n    seq.add(tf.keras.layers.ConvLSTM2D(16, (3, 3), padding=\"same\", return_sequences=True))\n    seq.add(tf.keras.layers.LayerNormalization())\n    # # # # #\n\n    seq.add(tf.keras.layers.TimeDistributed(\n        tf.keras.layers.Conv2DTranspose(16, (5, 5), strides=2, padding=\"same\")))\n    seq.add(tf.keras.layers.LayerNormalization())\n    seq.add(tf.keras.layers.TimeDistributed(\n        tf.keras.layers.Conv2DTranspose(32, (7, 7), strides=2, padding=\"same\")))\n    seq.add(tf.keras.layers.LayerNormalization())\n    seq.add(tf.keras.layers.TimeDistributed(\n        tf.keras.layers.Conv2DTranspose(64, (11, 11), strides=2, padding=\"same\")))\n    seq.add(tf.keras.layers.LayerNormalization())\n    seq.add(tf.keras.layers.TimeDistributed(\n        tf.keras.layers.Conv2D(1, (11, 11), activation=\"sigmoid\", padding=\"same\")))\n    \n    seq.compile(\n        loss='binary_crossentropy' if segment else 'mse', \n        optimizer=tf.keras.optimizers.Adam(lr=1e-4, decay=1e-5, epsilon=1e-6))\n    \n    return seq","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Also write a function to be able to fetch the tabular data we want the model to be able to predict:"},{"metadata":{"trusted":true},"cell_type":"code","source":"class DataGen(tf.keras.utils.Sequence):\n    def __init__(\n        self, patients, batch_size=4, unsupervised=True, scans=None, loc='train', \n        targets=['Recovering', 'slope'], augment=True, maxsize=10):\n        \n        '''\n        What augment parameter does?\n            Since we are using strided image slices, for each slice we could chose\n            the previous or the next slice randomly without causing much problem for \n            the model. This idea is exploited to create augmented scan data.\n            \n        We could use this for test data to maybe predict the uncertainity. But as I said,\n        unimplemented ideas :)\n        '''\n        \n        self.patients = patients\n        self.patient_count = len(patients)\n        self.batch_size = batch_size\n        self.unsupervised = unsupervised\n        self.loc = loc\n        self.targets = targets\n        self.scans = scans\n        self.maxsize = maxsize\n        self.augment = augment\n        \n    def __len__(self):\n        return self.patient_count // self.batch_size\n    \n    def __getitem__(self, idx):\n        \n        scan_slices, tab = [], []\n        for pat in self.patients[idx * self.batch_size: (idx+1) * self.batch_size]:\n            \n            if self.loc == 'train':\n                temp = np.load(scans[pat])\n            else:\n                temp = np.load(test_scans[pat])\n                \n            mask = np.linspace(0, len(temp)-1, self.maxsize).astype(int)\n                \n            if self.augment:\n                # possb combinations -> 3 ** (self.maxsize - 2)\n                mask = (mask + np.concatenate([[0], np.random.randint(-1, 2, size=self.maxsize-2), [0]]))\n                \n            scan_slices.append(temp[mask])\n                \n            if not self.unsupervised:\n                tab.append(self.fetch_tab(pat, targets=self.targets))\n        \n        scan_slices = np.stack(scan_slices)\n        \n        if self.unsupervised:\n            return scan_slices, scan_slices\n        else:\n            tab = np.stack(tab).astype(float)\n            return scan_slices, [tab[..., i] for i in range(len(self.targets))]\n        \n    def fetch_tab(self, pat, targets, data=base):\n        '''Given patient location, fetch the given as targets'''\n        return base.loc[base.Patient == pat, targets].values\n        \n    def __on_epoch_end__(self):\n        np.random.shuffle(self.patients)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"tf.keras.backend.clear_session()\n\nimgsize = 384\ntimesteps = 8\nbatch_size = 4\nsegment = True\ntrained = True   \n\nautoenc = get_model_small(imgsize=imgsize, timesteps=timesteps, segment=segment)\nprint (\"Features Per Patient: {:20}\\nTotal Number of Model Parameters: {:8}\"\n       .format(np.prod(autoenc.get_layer(\"enc_op\").output.shape[1:]), \n              autoenc.count_params()))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Loading the scan each and every time for the model training purpose will prove a bottle neck for model training. So the best solution would be to process them and save them as .npy files and load them as needed. \n\n`Warning:` Choosing higher imgsize or timesteps would cause an OSError resulting from insufficient physical memory space. They have already been optimized."},{"metadata":{"trusted":true},"cell_type":"code","source":"import os\nextracted_loc = \"./patient_saves\"\nif not os.path.exists(extracted_loc):\n    os.mkdir(extracted_loc)\n\nscans = {}\nfor pat in tqdm(train.Patient.unique()):\n    \n    temp = tf.image.resize(\n        extract_patient(f\"{main_dir}/train/{pat}\", max_size=timesteps*2)[..., tf.newaxis], \n        (imgsize, imgsize), method='nearest').numpy() / 2048\n\n    if segment:\n        temp = segment_scans(temp)\n\n    np.save(f\"{extracted_loc}/{pat}.npy\", temp)\n\n    scans[pat] = f\"{extracted_loc}/{pat}.npy\"\n\ntest_scans = {}\nfor pat in tqdm(test.Patient.unique()):\n    \n    temp = tf.image.resize(\n        extract_patient(f\"{main_dir}/test/{pat}\", max_size=timesteps)[..., tf.newaxis],\n        (imgsize, imgsize), method='nearest').numpy() / 2048\n\n    if segment:\n        temp = segment_scans(temp)\n\n    np.save(f\"{extracted_loc}/{pat}.npy\", temp)\n\n    test_scans[pat] = f\"{extracted_loc}/{pat}.npy\"\n\nprint (len(scans), len(test_scans))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Is our pipe Working good?"},{"metadata":{"trusted":true},"cell_type":"code","source":"%%capture\n\ntemp =  DataGen(train_patients, batch_size=batch_size, scans=scans, maxsize=timesteps)[0][0][0]\ntemp = tf.squeeze(temp)\nani = display_scan(temp)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print (temp.shape, tf.reduce_max(temp).numpy(), tf.reduce_min(temp).numpy())\nHTML(ani)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"if not trained:\n    hist = autoenc.fit(\n        DataGen(train_patients, batch_size=batch_size, scans=scans, maxsize=timesteps),\n        validation_data=DataGen(val_patients, batch_size=batch_size, \n                                scans=scans, maxsize=timesteps, augment=False),  \n        epochs=25, callbacks=[\n            tf.keras.callbacks.EarlyStopping(patience=3, mode='min', restore_best_weights=True),\n            tf.keras.callbacks.ModelCheckpoint(\n                'Model_Save-{val_loss:.3f}.hdf5', \n                save_best_only=True, mode='min')])\n    \n    # load the best weights\n    best_weight = min(glob.glob(\"./*.hdf5\"), key=lambda x: float(x.split(\"-\")[-1][:-5]))\n    \nelse:\n    best_weight = '../input/osic-weights/Model_Save-0.135.hdf5'\n    \n# load the best weights\nautoenc.load_weights(best_weight)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# validation embeddings loss\nautoenc.evaluate(DataGen(val_patients, batch_size=batch_size, scans=scans, maxsize=timesteps, augment=False))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"How well does our model make the predictions on unseen samples?"},{"metadata":{"trusted":true},"cell_type":"code","source":"%%capture\ntemp = np.load(scans[np.random.choice(val_patients)])\nani = display_scan(np.squeeze(autoenc.predict(\n    temp[np.linspace(0, len(temp)-1, timesteps).astype(int)][np.newaxis, :]\n)))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# the encoded part looks like this\nHTML(ani)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"*Yes, embeddings are pretty bad xD. Maybe could be improved by training for a longer epochs or using a more powerful model.*\n\nFor now let's use this embeddings to predict slopes for our patients. We create another model for this purpose. We feed in the embeddings and the model outputs the slope for that particular patient:"},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_slope_model(base_model, best_weight=best_weight):\n    \n    base_m = tf.keras.models.clone_model(base_model)\n    base_m.load_weights(best_weight)\n    \n    base_m.trainable = False\n    inp = base_m.input\n    enc_op = base_m.get_layer(\"enc_op\").output\n    enc_avg = tf.keras.layers.TimeDistributed(tf.keras.layers.GlobalAveragePooling2D())(enc_op)\n\n    lstm1 = tf.keras.layers.LSTM(\n        units=128, activation='relu',\n        kernel_initializer='he_normal')(enc_avg)\n\n    s_op = tf.keras.layers.Dense(1, name='Slope_op')(lstm1)\n    r_op = tf.keras.layers.Dense(1, activation='sigmoid', name='Recovery_op')(s_op)\n\n    model = tf.keras.models.Model(inputs=inp, outputs=[r_op, s_op])\n    model.compile(loss=['binary_crossentropy', 'mse'], optimizer='adam')\n    return model","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"folds = 3\ntotal_patients = train.Patient.unique()\nnp.random.shuffle(total_patients)\nval_len = len(total_patients) // folds\n\nval_preds, test_preds = np.zeros((val_len, 2)), np.zeros((len(test), 2))\n\nfor i in range(folds):\n    \n    model = get_slope_model(base_model=autoenc)\n   \n    val_patients = total_patients[(i)*val_len:(i+1)*val_len]\n    train_patients = np.setdiff1d(total_patients, val_patients)\n    \n    assert len(np.intersect1d(val_patients, train_patients)) == 0\n    \n    print (f\"Starting Fold #{i+1} {'='*36}>\")\n    \n    hist = model.fit(\n        DataGen(\n            train_patients, unsupervised=False, scans=scans, \n            batch_size=batch_size, maxsize=timesteps), \n        \n        validation_data=DataGen(\n            val_patients, unsupervised=False, scans=scans, \n            batch_size=batch_size, maxsize=timesteps, augment=False),\n        \n        verbose=0, epochs=50, callbacks=[tf.keras.callbacks.EarlyStopping(\n            patience=5, mode='min', monitor='val_loss', verbose=1, restore_best_weights=True)]\n    )\n    \n    print (\"Last Val Score (Pre Fine tune) : {:.2f}\".format(hist.history['val_loss'][-1]))\n\n    # fine tuning\n    model.trainable = True\n    model.compile(loss=['binary_crossentropy', 'mse'], optimizer=tf.keras.optimizers.Adam(0.0001))\n\n    hist = model.fit(\n        DataGen(train_patients, unsupervised=False, batch_size=batch_size, \n                scans=scans, maxsize=timesteps),\n        \n        validation_data=DataGen(val_patients, batch_size=batch_size, unsupervised=False, \n                                scans=scans, maxsize=timesteps, augment=False),\n        \n        verbose=0, initial_epoch=50, epochs=60,  callbacks=[tf.keras.callbacks.EarlyStopping(\n            patience=3, mode='min', monitor='val_loss', restore_best_weights=True, verbose=1)]\n    )\n    \n    print (\"Last Val Score (Post Fine tune): {:.2f}\\n{}\\n\".format(hist.history['val_loss'][-1], \"=\"*55))\n    \n    val_preds = val_preds + np.hstack(model.predict(DataGen(\n        val_patients, batch_size=1, scans=scans, maxsize=timesteps, augment=False))) / folds\n\n    test_preds = test_preds + np.hstack(model.predict(DataGen(\n        test['Patient'].tolist(), loc='test', batch_size=1, \n        scans=test_scans, maxsize=timesteps, augment=False))) / folds","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"THRESH = 0.5\ntest['Recovering'] = test_preds[:, 0] > THRESH\ntest['slope'] = test_preds[:, 1]\n\ntest.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"temp = pd.DataFrame(\n    np.stack([val_patients, val_preds[:, 1], val_preds[:, 0] > THRESH], axis=1),\n    columns=['Patient', 'slope', 'Recovering'])\n\ntemp = pd.concat([temp, base.loc[base.Patient.isin(train_patients), ['Patient', 'slope', 'Recovering']]])\ntemp[['slope', 'Recovering']] = temp[['slope', 'Recovering']].astype(float)\n\ntemp = train.merge(temp, on='Patient')\n\ntemp.loc[temp.Patient.isin(val_patients), 'Where'] = 'Val'\ntemp.loc[temp.Patient.isin(train_patients), 'Where'] = 'Train'\n\ntemp = base_shift(temp, q=0)\ntemp['factor'] = temp['Base_FVC'] / temp['Base_Percent']\ntemp['magic'] = temp['Base_FVC'] + (temp['Week_Offset'] * temp['slope'])\n\ntrain_cols = ['magic', 'Week_Offset', 'factor']\nlr = LinearRegression()\n\nlr.fit(\n    temp.loc[temp.Where == 'Train', train_cols], \n    temp.loc[temp.Where == 'Train', 'FVC']\n)\n\nprint (\"Train Score: {:.2f}\\nVal Score  : {:.2f}\".format(\n    laplace_log_likelihood(\n        temp.loc[temp.Where == 'Train', 'FVC'], \n        lr.predict(temp.loc[temp.Where == 'Train', train_cols]), 70),\n    \n    laplace_log_likelihood(\n    temp.loc[temp.Where == 'Val', 'FVC'], \n    lr.predict(temp.loc[temp.Where == 'Val', train_cols]), 70)    \n))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# creating x_test for predictions simaltanesly\nx_test = sub[['Patient', 'Weeks']].merge(\n    test.rename({\"Weeks\": \"Base_Week\", \n                 \"FVC\": \"Base_FVC\", \n                 \"Percent\": \"Base_Percent\"}, axis=1), \n    on='Patient')\n\n# create required features for predictions\nx_test['Week_Offset'] = x_test['Weeks'] - x_test['Base_Week']\nx_test['Percent'] = x_test['Base_Percent']\nx_test['factor'] = x_test['Base_FVC'] / x_test['Base_Percent']\nx_test['magic'] = x_test['Base_FVC'] + (x_test['slope'] * x_test['Week_Offset'])\n\nsub['FVC'] = lr.predict(x_test[train_cols])\nsub['Confidence'] = 200\n\n# save the submissions\nsub[['Patient_Week', 'FVC', 'Confidence']].to_csv(\"image_data_pred.csv\", index=False)\nsub.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Final Words:\n\n- I did try quantile regressions using Pinball loss. Didn't give me satisfactory results.\n- I tried my own neural network architectures and it performed poorly. In my own words: The problem is that it compromises on the RMSE instead of focusing on both, which is pretty apparent from the high Confidence values.\n\nIdeas left unimplemented:\n- The way doctors analyze scans is to compare the scans with similar scans that they had already seen (in their memory). So we could compare the scans using embeddings we had created earlier and for a patient, choose a patient most similar to and use the slope we had seen from him.\n- Using embeddings to measure similarity of patients is so much more reliable than using their slopes, pms_avg, etc. So this could be exploited for creating tabular data augmentations as well.\n\nThis concludes this Notebook. I hope you enjoyed reading this notebook of mine. Have a fantastic day, you are awesome <3"}],"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}