{"cells":[{"metadata":{},"cell_type":"markdown","source":"<a id=\"overview\"></a>\n# Overview 🧐\n[Pulmonary fibrosis is a lung disease that occurs when lung tissue becomes damaged and scarred. This thickened, stiff tissue makes it more difficult for your lungs to work properly. As pulmonary fibrosis worsens, you become progressively more short of breath.](https://www.mayoclinic.org/diseases-conditions/pulmonary-fibrosis/symptoms-causes/syc-20353690)<br>\n<font color=\"RoyalBlue\">肺線維症とは、肺の組織が傷つき、傷跡が残ることで起こる肺の病気です。この肥厚した硬い組織は、肺の正常な動作を困難にします。肺線維症が悪化すると、徐々に息切れがひどくなります。</font>\n<img src='https://i.imgur.com/edKPRik.png' width=\"600\">\nIn \"OSIC Pulmonary Fibrosis Progression\", we needs to predict a patient’s severity of decline in lung function based on a CT scan of their lungs by using AI machine learning. In detail, we must predict both a Forced vital capacity (FVC) and a confidence measure for each patient.<br>\n<font color=\"RoyalBlue\">「OSIC 肺線維症の進行」では、肺のCTスキャンに基づいて患者の肺機能の低下の重症度を AI 機械学習を用いて予測する必要があります。</font><br>\n<font color=\"RoyalBlue\">詳しく言うと、各患者の努力肺活量（FVC）と信頼度の両方を予測しなければなりません。</font>\n\n# Table of contents 📖\n* [Overview 🧐](#overview)\n* [Acknowledgements 🙇](#acknowledgements)\n* [Setup 💻](#setup)\n* [Load the data 📃](#load)\n* [Explore CSV data 📊](#explore)\n    * [Distribution of unique patients data 😷 (Age, Sex, SmokingStatus)](#unique)\n    * [Weeks distribution 📅](#weeks)\n    * [FVC & Percent distribution 💨](#percent)\n    * [Relationships between FVC and other variables 🤝](#fvc)\n* [Linear Decay (based on EfficientNets) 📷](#efficient)\n* [Multiple Quantile Regression 🌒](#quantile)\n    * [Data preprocessing for Multiple Quantile Regression 🧹](#quantile_d)\n    * [Build the model 🧠](#quantile_m)\n    * [Cross validation 💭](#quantile_c)\n* [Ensemble & Submit 📝](#submit)\n\n<a id=\"acknowledgements\"></a>\n# Acknowledgements 🙇\n\n- Ulrich GOUE's [Osic-Multiple-Quantile-Regression-Starter](https://www.kaggle.com/ulrich07/osic-multiple-quantile-regression-starter)\n- Michael Kazachok's [Linear Decay (based on ResNet CNN)](https://www.kaggle.com/miklgr500/linear-decay-based-on-resnet-cnn)\n- Wei Hao Khoong's [](http://)[EfficientNets + Quantile Regression (Inference)](https://www.kaggle.com/khoongweihao/efficientnets-quantile-regression-inference)\n\n<a id=\"setup\"></a>\n# Setup 💻\nAll seed values are fixed at 42.<br>\n<font color=\"RoyalBlue\">シード値は全て42で固定しています。</font><br>"},{"metadata":{"_kg_hide-output":true,"trusted":true},"cell_type":"code","source":"!pip install ../input/kerasapplications/keras-team-keras-applications-3b180cb -f ./ --no-index\n!pip install ../input/efficientnet/efficientnet-1.1.0/ -f ./ --no-index\n\nimport efficientnet.tfkeras as efn","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import os\nimport random\nimport pandas as pd\nimport numpy as np\nfrom tqdm.notebook import tqdm\nimport matplotlib.pyplot as plt\n%matplotlib inline\nimport seaborn as sns\n\nfrom sklearn.metrics import mean_absolute_error\nfrom sklearn.model_selection import KFold, train_test_split \nimport tensorflow as tf\nfrom tensorflow.keras import Model, backend\nimport tensorflow.keras.layers as L\nimport tensorflow.keras.models as M\nfrom tensorflow.keras.utils import Sequence\nfrom keras.utils.vis_utils import plot_model\n\nimport pydicom\nimport cv2\n\ndef seed_everything(seed=2020):\n    random.seed(seed)\n    os.environ['PYTHONHASHSEED'] = str(seed)\n    np.random.seed(seed)\n    tf.random.set_seed(seed)\n    \nseed_everything(42)\n\nconfig = tf.compat.v1.ConfigProto()\nconfig.gpu_options.allow_growth = True\nsession = tf.compat.v1.Session(config=config)\n\npalette_ro = [\"#ee2f35\", \"#fa7211\", \"#fbd600\", \"#75c731\", \"#1fb86e\", \"#0488cf\", \"#7b44ab\"]\n\nROOT = \"../input/osic-pulmonary-fibrosis-progression/\"","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"<a id=\"load\"></a>\n# Load the data 📃"},{"metadata":{"trusted":true},"cell_type":"code","source":"train = pd.read_csv(ROOT + \"train.csv\")\ntest = pd.read_csv(ROOT + \"test.csv\")\nsub = pd.read_csv(ROOT + \"sample_submission.csv\")\n\nprint(\"Training data shape: \", train.shape)\nprint(\"Test data shape: \", test.shape)\n\ntrain.head(10)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"* `Patient` - 各患者固有の ID（患者の DICOM フォルダの名前でもあります）\n* `Weeks` - ベースライン CT の前後の相対的な週数（負数の場合も）\n* `FVC` - 記録された努力肺活量（mL）\n* `Percent` - %FVC, パーセント肺活量。年齢・身長・性別から計算した予測 FVC に対する実際の FVC の割合\n* `Age` - 年齢\n* `Sex` - 性別（`Male` / `Female`）\n* `SmokingStatus` - 喫煙状態（`Never smoked` / `Ex-smoker` / `Currently smokes`）\n\n<a id=\"explore\"></a>\n# Explore CSV data 📊"},{"metadata":{"trusted":true},"cell_type":"code","source":"train.isnull().sum()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test.isnull().sum()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"There are no missing values in both `train` and `test`.<br>\n<font color=\"RoyalBlue\">train と test の両方に欠損値はありません。</font><br>\nBefore EDA, we will try to find duplicate rows in the `train` where the `Patient` and `Weeks` elements match.<br>\n<font color=\"RoyalBlue\">EDA の前に、train の中にある Patient と Weeks の要素が一致し重複している行を探してみます。</font>"},{"metadata":{"trusted":true},"cell_type":"code","source":"dupRows_train = train[train.duplicated(subset=['Patient', 'Weeks'], keep=False)]\n\nprint(\"There are {} duplicate rows here ({:.2f} percent of the total).\".format(len(dupRows_train), len(dupRows_train)/len(train)*100))\ndupRows_train","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"There don't seem to be too many of them. These duplicate rows should be removed.<br>\n<font color=\"RoyalBlue\">数はあまり多くないようです。これらの重複した行は削除しておきましょう。</font>"},{"metadata":{"trusted":true},"cell_type":"code","source":"train.drop_duplicates(subset=['Patient', 'Weeks'], keep=False, inplace=True)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The following table from [Is this Malware? [EDA, FE and lgb][updated]](https://www.kaggle.com/artgor/is-this-malware-eda-fe-and-lgb-updated)<br>\n\n<font color=\"RoyalBlue\">カラム名 / カラムごとのユニーク値数 / 最も出現頻度の高い値 / 最も出現頻度の高い値の出現回数 / 欠損損値の割合 / 最も多いカテゴリの割合 / dtypes を表示しています。<br>\ntrain における Patient の固有 ID 数は176のようです。</font>"},{"metadata":{"trusted":true},"cell_type":"code","source":"stats = []\nfor col in train.columns:\n    stats.append((col,\n                  train[col].nunique(),\n                  train[col].value_counts().index[0],\n                  train[col].value_counts().values[0],\n                  train[col].isnull().sum() * 100 / train.shape[0],\n                  train[col].value_counts(normalize=True, dropna=False).values[0] * 100,\n                  train[col].dtype))\nstats_df = pd.DataFrame(stats, columns=['Feature', 'Unique values', 'Most frequent item', 'Freuquence of most frequent item', 'Percentage of missing values', 'Percentage of values in the biggest category', 'Type'])\nstats_df.sort_values('Percentage of missing values', ascending=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"<a id=\"unique\"></a>\n## Distribution of unique patients data 😷 (Age, Sex, SmokingStatus)\n<font color=\"RoyalBlue\">では、train における Patient の固有 ID ごとの年齢、性別、喫煙状況の分布から見ていきましょう。</font>"},{"metadata":{"trusted":true},"cell_type":"code","source":"data = train.groupby(\"Patient\").first().reset_index(drop=True)\ndata.head()","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"fig, ((ax1, ax2, ax3), (ax4, ax5, ax6)) = plt.subplots(2, 3, figsize=(16, 12))\n\nsns.distplot(data[\"Age\"], ax=ax1, bins=data[\"Age\"].max()-data[\"Age\"].min()+1, color=palette_ro[1])\nax1.annotate(\"Min: {:,}\".format(data[\"Age\"].min()), xy=(data[\"Age\"].min(), 0.005), \n             xytext=(data[\"Age\"].min()-8, 0.02),\n             bbox=dict(boxstyle=\"round\", fc=\"none\", ec=\"gray\"),\n             arrowprops=dict(arrowstyle=\"->\",\n                             connectionstyle=\"arc3, rad=0.2\"))\nax1.annotate(\"Max: {:,}\".format(data[\"Age\"].max()), xy=(data[\"Age\"].max(), 0.005), \n             xytext=(data[\"Age\"].max()-2, 0.02),\n             bbox=dict(boxstyle=\"round\", fc=\"none\", ec=\"gray\"),\n             arrowprops=dict(arrowstyle=\"->\",\n                             connectionstyle=\"arc3, rad=-0.2\"))\nax1.axvline(x=data[\"Age\"].median(), color=palette_ro[0], linestyle=\"--\", alpha=0.5)\nax1.annotate(\"Med: {:.0f}\".format(data[\"Age\"].median()), xy=(data[\"Age\"].median(), 0.056), \n             xytext=(data[\"Age\"].median()-15, 0.065),\n             bbox=dict(boxstyle=\"round\", fc=\"none\", ec=\"gray\"),\n             arrowprops=dict(arrowstyle=\"->\",\n                             connectionstyle=\"arc3, rad=0.25\"))\n\nsns.countplot(x=\"Sex\", ax=ax2, data=data, palette=palette_ro[-2::-4])\nsns.countplot(x=\"SmokingStatus\", ax=ax3, data=data,\n              order=[\"Never smoked\", \"Ex-smoker\", \"Currently smokes\"], palette=palette_ro[-3::-2])\n\nsns.distplot(data[data[\"Sex\"]==\"Male\"].Age, label=\"Male\", ax=ax4, hist=False, color=palette_ro[5])\nsns.distplot(data[data[\"Sex\"]==\"Female\"].Age, label=\"Female\", ax=ax4, hist=False, color=palette_ro[1])\n\nsns.distplot(data[data[\"SmokingStatus\"]==\"Never smoked\"].Age, label=\"Never smoked\", ax=ax5, hist=False, color=palette_ro[4])\nsns.distplot(data[data[\"SmokingStatus\"]==\"Ex-smoker\"].Age, label=\"Ex-smoker\", ax=ax5, hist=False, color=palette_ro[2])\nsns.distplot(data[data[\"SmokingStatus\"]==\"Currently smokes\"].Age, label=\"Currently smokes\", ax=ax5, hist=False, color=palette_ro[0])\n\nsns.countplot(x=\"SmokingStatus\", ax=ax6, data=data, hue=\"Sex\",\n              order=[\"Never smoked\", \"Ex-smoker\", \"Currently smokes\"], palette=palette_ro[-2::-4])\n\nfig.suptitle(\"Distribution of unique patients data\", fontsize=18);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Compared to `male`, `female` seems to be of a wider age range and is less likely to smoke. And, `Never smoked` tend to be younger than `Ex-smoker`.<br>\n<font color=\"RoyalBlue\">男性に比べて女性は年齢層が幅広く、喫煙者が少ないようです。また、喫煙未経験者は元喫煙者よりも若い傾向にあります。</font>\n\n<a id=\"weeks\"></a>\n## Weeks distribution 📅"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(figsize=(16, 6))\n\nsns.distplot(train[\"Weeks\"], ax=ax, color=palette_ro[1], bins=train[\"Weeks\"].max()-train[\"Weeks\"].min()+1)\nax.annotate(\"Min: {:,}\".format(train[\"Weeks\"].min()), xy=(train[\"Weeks\"].min(), 0.005), \n            xytext=(train[\"Weeks\"].min()-8, 0.008),\n            bbox=dict(boxstyle=\"round\", fc=\"none\", ec=\"gray\"),\n            arrowprops=dict(arrowstyle=\"->\",\n                            connectionstyle=\"arc3, rad=0.2\"))\nax.annotate(\"Max: {:,}\".format(train[\"Weeks\"].max()), xy=(train[\"Weeks\"].max(), 0.005), \n            xytext=(train[\"Weeks\"].max()-2, 0.008),\n            bbox=dict(boxstyle=\"round\", fc=\"none\", ec=\"gray\"),\n            arrowprops=dict(arrowstyle=\"->\",\n                            connectionstyle=\"arc3, rad=-0.2\"))\nax.axvline(x=0, color=palette_ro[5], linestyle=\"--\", alpha=0.5)\nax.annotate(\"CT Scan\", xy=(0, 0.013), \n            xytext=(-12, 0.016),\n            bbox=dict(boxstyle=\"round\", fc=\"none\", ec=\"gray\"),\n            arrowprops=dict(arrowstyle=\"->\",\n                            connectionstyle=\"arc3, rad=0.2\"))\nax.axvline(x=train[\"Weeks\"].median(), color=palette_ro[0], linestyle=\"--\", alpha=0.5)\nax.annotate(\"Med: {:.0f}\".format(train[\"Weeks\"].median()), xy=(train[\"Weeks\"].median(), 0.020), \n            xytext=(train[\"Weeks\"].median()+2, 0.024),\n            bbox=dict(boxstyle=\"round\", fc=\"none\", ec=\"gray\"),\n            arrowprops=dict(arrowstyle=\"->\",\n                           connectionstyle=\"arc3, rad=-0.2\"))\n\nax.set_title(\"Weeks Distribution\", fontsize=18);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"<a id=\"percent\"></a>\n## FVC & Percent distribution 💨"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 12))\n\nsns.distplot(train[\"FVC\"], ax=ax1, color=palette_ro[5], hist=False)\nax1.annotate(\"Min: {:,}\".format(train[\"FVC\"].min()), xy=(train[\"FVC\"].min(), 0.00005), \n             xytext=(train[\"FVC\"].min()-300, 0.0001),\n             bbox=dict(boxstyle=\"round\", fc=\"none\", ec=\"gray\"),\n             arrowprops=dict(arrowstyle=\"->\",\n                            connectionstyle=\"arc3, rad=0.2\"))\nax1.annotate(\"Max: {:,}\".format(train[\"FVC\"].max()), xy=(train[\"FVC\"].max(), 0.00005), \n             xytext=(train[\"FVC\"].max()-200, 0.0001),\n             bbox=dict(boxstyle=\"round\", fc=\"none\", ec=\"gray\"),\n             arrowprops=dict(arrowstyle=\"->\",\n                             connectionstyle=\"arc3, rad=-0.2\"))\nax1.axvline(x=train[\"FVC\"].median(), color=palette_ro[0], linestyle=\"--\", alpha=0.5)\nax1.annotate(\"Med: {:,.0f}\".format(train[\"FVC\"].median()), xy=(train[\"FVC\"].median(), 0.00005), \n             xytext=(train[\"FVC\"].median()-750, 0.0001),\n             bbox=dict(boxstyle=\"round\", fc=\"none\", ec=\"gray\"),\n             arrowprops=dict(arrowstyle=\"->\",\n                           connectionstyle=\"arc3, rad=0.25\"))\n\nax1.set_title(\"FVC Distribution\", fontsize=16);\n\nsns.distplot(train[\"Percent\"], ax=ax2, color=palette_ro[3], hist=False)\nax2.annotate(\"Min: {:.2f}\".format(train[\"Percent\"].min()), xy=(train[\"Percent\"].min(), 0.0015), \n             xytext=(train[\"Percent\"].min()-8, 0.0040),\n             bbox=dict(boxstyle=\"round\", fc=\"none\", ec=\"gray\"),\n             arrowprops=dict(arrowstyle=\"->\",\n                            connectionstyle=\"arc3, rad=0.2\"))\nax2.annotate(\"Max: {:.2f}\".format(train[\"Percent\"].max()), xy=(train[\"Percent\"].max(), 0.0015), \n             xytext=(train[\"Percent\"].max()-4, 0.0040),\n             bbox=dict(boxstyle=\"round\", fc=\"none\", ec=\"gray\"),\n             arrowprops=dict(arrowstyle=\"->\",\n                             connectionstyle=\"arc3, rad=-0.2\"))\nax2.axvline(x=train[\"Percent\"].median(), color=palette_ro[0], linestyle=\"--\", alpha=0.5)\nax2.annotate(\"Med: {:.2f}\".format(train[\"Percent\"].median()), xy=(train[\"Percent\"].median(), 0.0015), \n             xytext=(train[\"Percent\"].median()-17, 0.0040),\n             bbox=dict(boxstyle=\"round\", fc=\"none\", ec=\"gray\"),\n             arrowprops=dict(arrowstyle=\"->\",\n                           connectionstyle=\"arc3, rad=0.25\"))\n\nax2.set_title(\"Percent Distribution\", fontsize=16);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"<a id=\"fvc\"></a>\n## Relationships between FVC and other variables 🤝\n\nLet's look at the relationships between the objective variable, `FVC`, and other variables.<br>\n<font color=\"RoyalBlue\">目的変数である FVC と他の変数との関係を見ていきましょう。</font>"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 12))\n\nsns.distplot(train[train[\"Sex\"]==\"Male\"].FVC, label=\"Male\", ax=ax1, hist=False, color=palette_ro[5])\nax1.axvline(x=train[train[\"Sex\"]==\"Male\"].FVC.median(), color=palette_ro[5], linestyle=\"--\", alpha=0.5)\nax1.annotate(\"Male\\nMed: {:,.0f}\".format(train[train[\"Sex\"]==\"Male\"].FVC.median()), xy=(train[train[\"Sex\"]==\"Male\"].FVC.median(), 0.0006), \n             xytext=(train[train[\"Sex\"]==\"Male\"].FVC.median()+100, 0.00065),\n             bbox=dict(boxstyle=\"round\", fc=\"none\", ec=\"gray\"),\n             arrowprops=dict(arrowstyle=\"->\",\n                           connectionstyle=\"arc3, rad=-0.2\"))\nsns.distplot(train[train[\"Sex\"]==\"Female\"].FVC, label=\"Female\", ax=ax1, hist=False, color=palette_ro[1])\nax1.axvline(x=train[train[\"Sex\"]==\"Female\"].FVC.median(), color=palette_ro[1], linestyle=\"--\", alpha=0.5)\nax1.annotate(\"Female\\nMed: {:,.0f}\".format(train[train[\"Sex\"]==\"Female\"].FVC.median()), xy=(train[train[\"Sex\"]==\"Female\"].FVC.median(), 0.0008), \n             xytext=(train[train[\"Sex\"]==\"Female\"].FVC.median()+100, 0.00085),\n             bbox=dict(boxstyle=\"round\", fc=\"none\", ec=\"gray\"),\n             arrowprops=dict(arrowstyle=\"->\",\n                           connectionstyle=\"arc3, rad=-0.2\"))\n\nsns.distplot(train[train[\"SmokingStatus\"]==\"Never smoked\"].FVC, label=\"Never smoked\", ax=ax2, hist=False, color=palette_ro[4])\nax2.axvline(x=train[train[\"SmokingStatus\"]==\"Never smoked\"].FVC.median(), color=palette_ro[4], linestyle=\"--\", alpha=0.5)\nax2.annotate(\"Never smoked\\nMed: {:.0f}\".format(train[train[\"SmokingStatus\"]==\"Never smoked\"].FVC.median()), xy=(train[train[\"SmokingStatus\"]==\"Never smoked\"].FVC.median(), 0.0005), \n             xytext=(train[train[\"SmokingStatus\"]==\"Never smoked\"].FVC.median()-1000, 0.00055),\n             bbox=dict(boxstyle=\"round\", fc=\"none\", ec=\"gray\"),\n             arrowprops=dict(arrowstyle=\"->\",\n                           connectionstyle=\"arc3, rad=0.25\"))\nsns.distplot(train[train[\"SmokingStatus\"]==\"Ex-smoker\"].FVC, label=\"Ex-smoker\", ax=ax2, hist=False, color=palette_ro[2])\nax2.axvline(x=train[train[\"SmokingStatus\"]==\"Ex-smoker\"].FVC.median(), color=palette_ro[2], linestyle=\"--\", alpha=0.75)\nax2.annotate(\"Ex-smoker\\nMed: {:.0f}\".format(train[train[\"SmokingStatus\"]==\"Ex-smoker\"].FVC.median()), xy=(train[train[\"SmokingStatus\"]==\"Ex-smoker\"].FVC.median(), 0.00058), \n             xytext=(train[train[\"SmokingStatus\"]==\"Ex-smoker\"].FVC.median()-1200, 0.0007),\n             bbox=dict(boxstyle=\"round\", fc=\"none\", ec=\"gray\"),\n             arrowprops=dict(arrowstyle=\"->\",\n                           connectionstyle=\"arc3, rad=-0.25\"))\nsns.distplot(train[train[\"SmokingStatus\"]==\"Currently smokes\"].FVC, label=\"Currently smokes\", ax=ax2, hist=False, color=palette_ro[0])\nax2.axvline(x=train[train[\"SmokingStatus\"]==\"Currently smokes\"].FVC.median(), color=palette_ro[0], linestyle=\"--\", alpha=0.5)\nax2.annotate(\"Currently smokes\\nMed: {:.0f}\".format(train[train[\"SmokingStatus\"]==\"Currently smokes\"].FVC.median()), xy=(train[train[\"SmokingStatus\"]==\"Currently smokes\"].FVC.median(), 0.0009), \n             xytext=(train[train[\"SmokingStatus\"]==\"Currently smokes\"].FVC.median()+400, 0.00095),\n             bbox=dict(boxstyle=\"round\", fc=\"none\", ec=\"gray\"),\n             arrowprops=dict(arrowstyle=\"->\",\n                           connectionstyle=\"arc3, rad=0.25\"))\n\nax1.set_title(\"Relationship between FVC and Sex\", fontsize=16)\nax2.set_title(\"Relationship between FVC and SmokingStatus\", fontsize=16);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"You can see that `FVC` of `Female` tends to be much lower than of `Male`, and there is also a difference in FVC by `SmokingStatus`, but this may be because there are more `Female` in `Never smoked`. Let's check it out.<br>\n<font color=\"RoyalBlue\">男性と比べると女性の FVC はかなり低くなる傾向にあることが分かります。喫煙状態によっても FVC に差が出ていますが、これは喫煙未経験者に女性が多いためかもしれません。実際に確認してみましょう。</font>"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"train_m = train[train[\"Sex\"]==\"Male\"].reset_index(drop=True)\ntrain_f = train[train[\"Sex\"]==\"Female\"].reset_index(drop=True)\n\nfig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 12))\n\nsns.distplot(train_m[train_m[\"SmokingStatus\"]==\"Never smoked\"].FVC, label=\"Never smoked\", ax=ax1, hist=False, color=palette_ro[4])\nax1.axvline(x=train_m[train_m[\"SmokingStatus\"]==\"Never smoked\"].FVC.median(), color=palette_ro[4], linestyle=\"--\", alpha=0.5)\nax1.annotate(\"Never smoked\\nMed: {:.0f}\".format(train_m[train_m[\"SmokingStatus\"]==\"Never smoked\"].FVC.median()), xy=(train_m[train_m[\"SmokingStatus\"]==\"Never smoked\"].FVC.median(), 0.0005), \n             xytext=(train_m[train_m[\"SmokingStatus\"]==\"Never smoked\"].FVC.median()-1400, 0.0006),\n             bbox=dict(boxstyle=\"round\", fc=\"none\", ec=\"gray\"),\n             arrowprops=dict(arrowstyle=\"->\",\n                           connectionstyle=\"arc3, rad=0.25\"))\nsns.distplot(train_m[train_m[\"SmokingStatus\"]==\"Ex-smoker\"].FVC, label=\"Ex-smoker\", ax=ax1, hist=False, color=palette_ro[2])\nax1.axvline(x=train_m[train_m[\"SmokingStatus\"]==\"Ex-smoker\"].FVC.median(), color=palette_ro[2], linestyle=\"--\", alpha=0.75)\nax1.annotate(\"Ex-smoker\\nMed: {:.0f}\".format(train_m[train_m[\"SmokingStatus\"]==\"Ex-smoker\"].FVC.median()), xy=(train_m[train_m[\"SmokingStatus\"]==\"Ex-smoker\"].FVC.median(), 0.00063), \n             xytext=(train_m[train_m[\"SmokingStatus\"]==\"Ex-smoker\"].FVC.median()-1400, 0.00045),\n             bbox=dict(boxstyle=\"round\", fc=\"none\", ec=\"gray\"),\n             arrowprops=dict(arrowstyle=\"->\",\n                           connectionstyle=\"arc3, rad=-0.2\"))\nsns.distplot(train_m[train_m[\"SmokingStatus\"]==\"Currently smokes\"].FVC, label=\"Currently smokes\", ax=ax1, hist=False, color=palette_ro[0])\nax1.axvline(x=train_m[train_m[\"SmokingStatus\"]==\"Currently smokes\"].FVC.median(), color=palette_ro[0], linestyle=\"--\", alpha=0.5)\nax1.annotate(\"Currently smokes\\nMed: {:.0f}\".format(train_m[train_m[\"SmokingStatus\"]==\"Currently smokes\"].FVC.median()), xy=(train_m[train_m[\"SmokingStatus\"]==\"Currently smokes\"].FVC.median(), 0.00066), \n             xytext=(train_m[train_m[\"SmokingStatus\"]==\"Currently smokes\"].FVC.median()+400, 0.00055),\n             bbox=dict(boxstyle=\"round\", fc=\"none\", ec=\"gray\"),\n             arrowprops=dict(arrowstyle=\"->\",\n                           connectionstyle=\"arc3, rad=0.25\"))\n\nsns.distplot(train_f[train_f[\"SmokingStatus\"]==\"Never smoked\"].FVC, label=\"Never smoked\", ax=ax2, hist=False, color=palette_ro[4])\nax2.axvline(x=train_f[train_f[\"SmokingStatus\"]==\"Never smoked\"].FVC.median(), color=palette_ro[4], linestyle=\"--\", alpha=0.5)\nax2.annotate(\"Never smoked\\nMed: {:.0f}\".format(train_f[train_f[\"SmokingStatus\"]==\"Never smoked\"].FVC.median()), xy=(train_f[train_f[\"SmokingStatus\"]==\"Never smoked\"].FVC.median(), 0.001), \n             xytext=(train_f[train_f[\"SmokingStatus\"]==\"Never smoked\"].FVC.median()-600, 0.0015),\n             bbox=dict(boxstyle=\"round\", fc=\"none\", ec=\"gray\"),\n             arrowprops=dict(arrowstyle=\"->\",\n                           connectionstyle=\"arc3, rad=-0.2\"))\nsns.distplot(train_f[train_f[\"SmokingStatus\"]==\"Ex-smoker\"].FVC, label=\"Ex-smoker\", ax=ax2, hist=False, color=palette_ro[2])\nax2.axvline(x=train_f[train_f[\"SmokingStatus\"]==\"Ex-smoker\"].FVC.median(), color=palette_ro[2], linestyle=\"--\", alpha=0.75)\nax2.annotate(\"Ex-smoker\\nMed: {:.0f}\".format(train_f[train_f[\"SmokingStatus\"]==\"Ex-smoker\"].FVC.median()), xy=(train_f[train_f[\"SmokingStatus\"]==\"Ex-smoker\"].FVC.median(), 0.0013), \n             xytext=(train_f[train_f[\"SmokingStatus\"]==\"Ex-smoker\"].FVC.median()+100, 0.0018),\n             bbox=dict(boxstyle=\"round\", fc=\"none\", ec=\"gray\"),\n             arrowprops=dict(arrowstyle=\"->\",\n                           connectionstyle=\"arc3, rad=-0.25\"))\nsns.distplot(train_f[train_f[\"SmokingStatus\"]==\"Currently smokes\"].FVC, label=\"Currently smokes\", ax=ax2, hist=False, color=palette_ro[0])\nax2.axvline(x=train_f[train_f[\"SmokingStatus\"]==\"Currently smokes\"].FVC.median(), color=palette_ro[0], linestyle=\"--\", alpha=0.5)\nax2.annotate(\"Currently smokes\\nMed: {:.0f}\".format(train_f[train_f[\"SmokingStatus\"]==\"Currently smokes\"].FVC.median()), xy=(train_f[train_f[\"SmokingStatus\"]==\"Currently smokes\"].FVC.median(), 0.0035), \n             xytext=(train_f[train_f[\"SmokingStatus\"]==\"Currently smokes\"].FVC.median()+200, 0.004),\n             bbox=dict(boxstyle=\"round\", fc=\"none\", ec=\"gray\"),\n             arrowprops=dict(arrowstyle=\"->\",\n                           connectionstyle=\"arc3, rad=0.25\"))\n\nax1.set_title(\"Relationship between FVC and SmokingStatus in Male\", fontsize=16)\nax2.set_title(\"Relationship between FVC and SmokingStatus in Female\", fontsize=16);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"When limited to `Male`, `FVC` does not seem to change much with `SmokingStatus`. In the case of `Female`, there is a difference, but this is probably due to the small sample size (especially for `Currently smokes`). It may also be important to consider that patients who are `Currently smokes` are likely to be less severely affected.<br>\n<font color=\"RoyalBlue\">男性に限定してみると、FVC は喫煙状態によってはあまり変化しないようです。女性の場合は差が出ていますが、これはサンプル数が少ないためでしょう（特に現喫煙者）。現喫煙者の患者には重症者が少ないであろうことも考慮すべきかもしれません。</font>"},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 12))\n\nsns.scatterplot(x=train[\"FVC\"], y=train[\"Age\"], ax=ax1,\n                palette=[palette_ro[5], palette_ro[1]], hue=train[\"Sex\"], style=train[\"Sex\"])\nsns.scatterplot(x=train[\"FVC\"], y=train[\"Age\"], ax=ax2,\n                palette=[palette_ro[2], palette_ro[4], palette_ro[0]], hue=train[\"SmokingStatus\"], style=train[\"SmokingStatus\"])\n\nfig.suptitle(\"Correlation between FVC and Age (Pearson Corr: {:.4f})\".format(train[\"FVC\"].corr(train[\"Age\"])), fontsize=16);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"There was almost no correlation between `FVC` and `Age`.<br>\n<font color=\"RoyalBlue\">FVC と年齢には相関はほぼ見られませんでした。</font>"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 12))\n\nsns.scatterplot(x=train[\"FVC\"], y=train[\"Weeks\"], ax=ax1,\n                palette=[palette_ro[5], palette_ro[1]], hue=train[\"Sex\"], style=train[\"Sex\"])\nsns.scatterplot(x=train[\"FVC\"], y=train[\"Weeks\"], ax=ax2,\n                palette=[palette_ro[2], palette_ro[4], palette_ro[0]], hue=train[\"SmokingStatus\"], style=train[\"SmokingStatus\"])\n\nfig.suptitle(\"Correlation between FVC and Weeks (Pearson Corr: {:.4f})\".format(train[\"FVC\"].corr(train[\"Weeks\"])), fontsize=16);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"There was almost no correlation between `FVC` and `Weeks` either.<br>\n<font color=\"RoyalBlue\">FVC と週数にも相関はほぼ見られませんでした。</font>"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 12))\n\nsns.scatterplot(x=train[\"FVC\"], y=train[\"Age\"], ax=ax1,\n                palette=[palette_ro[5], palette_ro[1]], hue=train[\"Sex\"], style=train[\"Sex\"])\nsns.scatterplot(x=train[\"FVC\"], y=train[\"Age\"], ax=ax2,\n                palette=[palette_ro[2], palette_ro[4], palette_ro[0]], hue=train[\"SmokingStatus\"], style=train[\"SmokingStatus\"])\n\nfig.suptitle(\"Correlation between FVC and Age (Pearson Corr: {:.4f})\".format(train[\"FVC\"].corr(train[\"Age\"])), fontsize=16);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"There was a positive correlation between `FVC` and `Percent`, which is not surprising since `Percent` is a value calculated from `FVC` and other data.<br>\n<font color=\"RoyalBlue\">FVC と Percent には正の相関が見られました。Percent は FVC などから算出する値なので当然といえば当然です。</font>\n\n<a id=\"efficient\"></a>\n# Linear Decay (based on EfficientNets) 📷\nFirst of all, let's try to make predictions from DICOM and other data. We make a function to put `Age`, `Sex`, and `SmokingStatus` into NumPy array.<br>\n<font color=\"RoyalBlue\">それでは、まずは DICOM を含めたデータから予測を出してみましょう。Age, Sex, SmokingStatus を変換しつつ NumPy array にまとめる関数を作成します。</font>"},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_tab(df):\n    vector = [(df.Age.values[0] - 30) / 30]\n    \n    if df.Sex.values[0] == 'male':\n       vector.append(0)\n    else:\n       vector.append(1)\n    \n    if df.SmokingStatus.values[0] == 'Never smoked':\n        vector.extend([0,0])\n    elif df.SmokingStatus.values[0] == 'Ex-smoker':\n        vector.extend([1,1])\n    elif df.SmokingStatus.values[0] == 'Currently smokes':\n        vector.extend([0,1])\n    else:\n        vector.extend([1,0])\n    return np.array(vector)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Next, we organize the data by Patient's unique ID.<br>\nIn the for statement, we extract `FVC` from `fvc` and `Weeks` from `weeks` in a NumPy array format, respectively.\nWe join weeks and the array with all 1's in the vertical direction and transpose them, and then assign them to `c`.\nThen find the least-squares of `c` and `fvc`, and add the slope to `A`.\nWe then add the value obtained from the above function `get_tab` to `TAB` and the unique ID of Patient to `p`.<br>\n<font color=\"RoyalBlue\">次に、Patient の固有 ID ごとにデータを整理します。<br>\nfor 文の中で fvc に FVC, weeks に Weeks をそれぞれ NumPy array 形式で取り出します。。<br>\nweeks と要素が全て1の配列を縦方向に結合し、転置したものを c とします。<br>\nそして c と fvc の最小二乗解を求め、その傾きを A に追加します。<br>\nさらに上記の上記の関数 get_tab で取得した値を TAB に、Patient の固有 ID を p に追加していきます。</font>"},{"metadata":{"trusted":true},"cell_type":"code","source":"A = {} \nTAB = {} \nP = [] \nfor i, p in tqdm(enumerate(train.Patient.unique())):\n    sub = train.loc[train.Patient == p, :] \n    fvc = sub.FVC.values\n    weeks = sub.Weeks.values\n    c = np.vstack([weeks, np.ones(len(weeks))]).T\n    a, b = np.linalg.lstsq(c, fvc)[0]\n    \n    A[p] = a\n    TAB[p] = get_tab(sub)\n    P.append(p)\n\nfig, ax = plt.subplots(1, 1, figsize=(16, 6))\nsns.distplot(list(A.values()), ax=ax, color=palette_ro[1]);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Create a function to read a DICOM file from the given path.\nFor ease of use in deep learning we will divide the value by 2048 and then crop and resize the image.<br>\n<font color=\"RoyalBlue\">渡されたパスから DICOM ファイルを読み込む関数を作成します。<br>\nディープラーニングで扱いやすくするために値を 2**11 で割り、さらに画像をクロップ・リサイズします。</font>"},{"metadata":{"trusted":true},"cell_type":"code","source":"# def get_img(path):\n#     d = pydicom.dcmread(path)\n#     return cv2.resize(d.pixel_array / 2**11, (528, 528))    # changed from 512\n\n# https://www.kaggle.com/allunia/pulmonary-dicom-preprocessing\ndef get_img(path, new_shape=(528, 528)):\n    d = pydicom.dcmread(path)\n    scan = d.pixel_array / 2**11\n    \n    left = int((scan.shape[0]-512)/2)\n    right = int((scan.shape[0]+512)/2)\n    top = int((scan.shape[1]-512)/2)\n    bottom = int((scan.shape[1]+512)/2)\n    \n    img = scan[top:bottom, left:right]\n    cropped_resized_scan = cv2.resize(img, new_shape, interpolation=cv2.INTER_LANCZOS4)\n    return cropped_resized_scan\n\n# get_img(\"../input/osic-pulmonary-fibrosis-progression/train/ID00007637202177411956430/1.dcm\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class IGenerator(Sequence):\n    BAD_ID = ['ID00011637202177653955184', 'ID00052637202186188008618']\n    def __init__(self, keys, a, tab, batch_size=32):\n        self.keys = [k for k in keys if k not in self.BAD_ID]\n        self.a = a\n        self.tab = tab\n        self.batch_size = batch_size\n        \n        self.train_data = {}\n        for p in train.Patient.values:\n            self.train_data[p] = os.listdir(f'../input/osic-pulmonary-fibrosis-progression/train/{p}/')\n    \n    def __len__(self):\n        return 1000\n    \n    def __getitem__(self, idx):\n        x = []\n        a, tab = [], [] \n        keys = np.random.choice(self.keys, size = self.batch_size)\n        for k in keys:\n            try:\n                i = np.random.choice(self.train_data[k], size=1)[0]\n                img = get_img(f'../input/osic-pulmonary-fibrosis-progression/train/{k}/{i}')\n                x.append(img)\n                a.append(self.a[k])\n                tab.append(self.tab[k])\n            except:\n                print(k, i)\n       \n        x,a,tab = np.array(x), np.array(a), np.array(tab)\n        x = np.expand_dims(x, axis=-1)\n        return [x, tab] , a","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's create a model.\nWe take image data and run the tensor through `EfficientNetB6` and `GlobalAveragePooling2D`, and we take the pre-processed CSV data and add Gaussian noise to the tensor and concatenate them together to output the model.\nThe model weights are trained.<br>\n<font color=\"RoyalBlue\">では、モデルを作成していきます。<br>\n画像データを受け取って EfficientNetB6 と GlobalAveragePooling2D に通したテンソルと、前処理済みの CSV データを受け取ってガウシアンノイズを加えたテンソルを連結して出力します。<br>\nモデルの重みは訓練済みのものを使用します。</font>"},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\ndef get_efficientnet(model, shape):\n    models_dict = {\n        'b0': efn.EfficientNetB0(input_shape=shape,weights=None,include_top=False),\n        'b1': efn.EfficientNetB1(input_shape=shape,weights=None,include_top=False),\n        'b2': efn.EfficientNetB2(input_shape=shape,weights=None,include_top=False),\n        'b3': efn.EfficientNetB3(input_shape=shape,weights=None,include_top=False),\n        'b4': efn.EfficientNetB4(input_shape=shape,weights=None,include_top=False),\n        'b5': efn.EfficientNetB5(input_shape=shape,weights=None,include_top=False),\n        'b6': efn.EfficientNetB6(input_shape=shape,weights=None,include_top=False),\n        'b7': efn.EfficientNetB7(input_shape=shape,weights=None,include_top=False)\n    }\n    return models_dict[model]\n\ndef build_model(shape=(528, 528, 1), model_class=None):    # changed from 512\n    inp = L.Input(shape=shape)\n    base = get_efficientnet(model_class, shape)\n    x = base(inp)\n    x = L.GlobalAveragePooling2D()(x)\n    inp2 = L.Input(shape=(4,))\n    x2 = L.GaussianNoise(0.2)(inp2)\n    x = L.Concatenate()([x, x2]) \n    x = L.Dropout(0.32)(x)    # changed from 0.4\n    x = L.Dense(1)(x)\n    model = Model([inp, inp2] , x)\n    \n    weights = [w for w in os.listdir('../input/osic-model-weights') if model_class in w][0]\n    model.load_weights('../input/osic-model-weights/' + weights)\n    return model\n\nmodel_classes = [\"b6\"] #['b0','b1','b2','b3',b4','b5','b6','b7']    # changed from b5\nmodels = [build_model(model_class=m, shape=(528, 528, 1)) for m in model_classes]    # changed from 512\nprint('Number of models: ' + str(len(models)))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_model(models[0])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We will create a modified version of the Laplace Log Likelihood function, which is the evaluation function for this competition.<br>\n<font color=\"RoyalBlue\">このコンペでの評価関数である修正版ラプラス対数尤度の関数を作成しておきます。</font>"},{"metadata":{"trusted":true},"cell_type":"code","source":"def score(fvc_true, fvc_pred, sigma):\n    sigma_clip = np.maximum(sigma, 70)\n    delta = np.abs(fvc_true - fvc_pred)\n    delta = np.minimum(delta, 1000)\n    sq2 = np.sqrt(2)\n    metric = (delta / sigma_clip)*sq2 + np.log(sigma_clip * sq2)\n    return np.mean(metric)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now, we will use the model to make predictions.<br>\n<font color=\"RoyalBlue\">それでは、モデルを使って予測を行います。</font>"},{"metadata":{"trusted":true},"cell_type":"code","source":"tr_p, vl_p = train_test_split(P, \n                              shuffle=True, \n                              train_size=0.8)\n\nsubs = []\nfor model in models:\n    metric = []\n    for q in tqdm(range(1, 10)):\n        m = []\n        for p in vl_p:\n            x = [] \n            tab = [] \n\n            if p in ['ID00011637202177653955184', 'ID00052637202186188008618']:\n                continue\n\n            ldir = os.listdir(f'../input/osic-pulmonary-fibrosis-progression/train/{p}/')\n            for i in ldir:\n                if int(i[:-4]) / len(ldir) < 0.8 and int(i[:-4]) / len(ldir) > 0.15:\n                    x.append(get_img(f'../input/osic-pulmonary-fibrosis-progression/train/{p}/{i}')) \n                    tab.append(get_tab(train.loc[train.Patient == p, :])) \n            if len(x) < 1:\n                continue\n            tab = np.array(tab) \n\n            x = np.expand_dims(x, axis=-1) \n            _a = model.predict([x, tab]) \n            a = np.quantile(_a, q / 10)\n\n            percent_true = train.Percent.values[train.Patient == p]\n            fvc_true = train.FVC.values[train.Patient == p]\n            weeks_true = train.Weeks.values[train.Patient == p]\n\n            fvc = a * (weeks_true - weeks_true[0]) + fvc_true[0]\n            percent = percent_true[0] - a * abs(weeks_true - weeks_true[0])\n            m.append(score(fvc_true, fvc, percent))\n        print(np.mean(m))\n        metric.append(np.mean(m))\n\n    q = (np.argmin(metric) + 1)/ 10\n\n    sub = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/sample_submission.csv') \n    test = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/test.csv') \n    A_test, B_test, P_test, W, FVC = {}, {}, {}, {}, {} \n    STD, WEEK = {}, {} \n    for p in test.Patient.unique():\n        x = [] \n        tab = [] \n        ldir = os.listdir(f'../input/osic-pulmonary-fibrosis-progression/test/{p}/')\n        for i in ldir:\n            if int(i[:-4]) / len(ldir) < 0.8 and int(i[:-4]) / len(ldir) > 0.15:\n                x.append(get_img(f'../input/osic-pulmonary-fibrosis-progression/test/{p}/{i}')) \n                tab.append(get_tab(test.loc[test.Patient == p, :])) \n        if len(x) <= 1:\n            continue\n        tab = np.array(tab) \n\n        x = np.expand_dims(x, axis=-1) \n        _a = model.predict([x, tab]) \n        a = np.quantile(_a, q)\n        A_test[p] = a\n        B_test[p] = test.FVC.values[test.Patient == p] - a*test.Weeks.values[test.Patient == p]\n        P_test[p] = test.Percent.values[test.Patient == p] \n        WEEK[p] = test.Weeks.values[test.Patient == p]\n\n    for k in sub.Patient_Week.values:\n        p, w = k.split('_')\n        w = int(w) \n\n        fvc = A_test[p] * w + B_test[p]\n        sub.loc[sub.Patient_Week == k, \"FVC\"] = fvc\n        sub.loc[sub.Patient_Week == k, \"Confidence\"] = (\n            P_test[p] - A_test[p] * abs(WEEK[p] - w) \n    ) \n\n    _sub = sub[[\"Patient_Week\",\"FVC\",\"Confidence\"]].copy()\n    subs.append(_sub)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"N = len(subs)\nsub = subs[0].copy() # ref\nsub[\"FVC\"] = 0\nsub[\"Confidence\"] = 0\nfor i in range(N):\n    sub[\"FVC\"] += subs[0][\"FVC\"] * (1/N)\n    sub[\"Confidence\"] += subs[0][\"Confidence\"] * (1/N)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sub.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sub[[\"Patient_Week\", \"FVC\", \"Confidence\"]].to_csv(\"submission_img.csv\", index=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"img_sub = sub[[\"Patient_Week\",\"FVC\",\"Confidence\"]].copy()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Linear Decay predictions have been made!<br>\n<font color=\"RoyalBlue\">Linear Decay での予測が作成できました！</font>\n\n<a id=\"quantile\"></a>\n# Multiple Quantile Regression 🌒\n<a id=\"quantile_d\"></a>\n## Data preprocessing for Multiple Quantile Regression 🧹\n\nNext, let's pre-process the data for multiple quantile regression. First, check the format of the `sub` (sample_submission.csv).<br>\n<font color=\"RoyalBlue\">次に、重分位点回帰のためのデータの前処理を行っていきましょう。まず、sub (sample_submission.csv) の形式をチェックします。</font>"},{"metadata":{"trusted":true},"cell_type":"code","source":"sub = pd.read_csv(ROOT + \"sample_submission.csv\")\nsub.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Split `Patient_Week` in `sub` into `Patient` and `Weeks`, according to the `train` and `test` formats. Then, attach the `Patient` to `Patient` in `sub` and merge it with the `Patient`. This makes it easier to handle the prediction.<br>\n<font color=\"RoyalBlue\">sub の Patient_Week を Patient と Weeks に分割し、train や test の形式に合わせます。そして、sub に test を Patient に紐づけて結合します。こうすれば、予測時に簡単に処理を行えるようになります。</font>"},{"metadata":{"trusted":true},"cell_type":"code","source":"sub['Patient'] = sub['Patient_Week'].apply(lambda x:x.split('_')[0])\nsub['Weeks'] = sub['Patient_Week'].apply(lambda x: int(x.split('_')[-1]))\nsub =  sub[['Patient', 'Weeks', 'Confidence', 'Patient_Week']]\nsub = sub.merge(test.drop('Weeks', axis=1), on=\"Patient\")\nsub.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Add `Where` column to all the dataframes.<br>\n<font color=\"RoyalBlue\">すべてのデータフレームに Where 列を追加します。</font><br>\nThen, in order to process the `train`, `test` and `sub` at the same time, these three are concatenated vertically into `data`.<br>\n<font color=\"RoyalBlue\">そして、train, test, sub を同時にデータ処理するために、縦方向に結合して data としておきます。</font>"},{"metadata":{"trusted":true},"cell_type":"code","source":"train['WHERE'] = 'train'\ntest['WHERE'] = 'val'\nsub['WHERE'] = 'test'\ndata = train.append([test, sub])\n\nprint(train.shape, test.shape, sub.shape, data.shape)\nprint(train.Patient.nunique(), test.Patient.nunique(), sub.Patient.nunique(), data.Patient.nunique())\n\ndata.head(10)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Add a `min_week` column for the minimum number of weeks per Patient.<br>\n<font color=\"RoyalBlue\">Patient ごとの最小の週数を示す min_week 列を追加します。</font>"},{"metadata":{"trusted":true},"cell_type":"code","source":"data['min_week'] = data['Weeks']\ndata.loc[data.WHERE=='test', 'min_week'] = np.nan\ndata['min_week'] = data.groupby('Patient')['min_week'].transform('min')\n\ndata.head(10)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"From here, calculate `base_FVC` (= the `FVC` of `Patient` at `min_week`) and the `base_week` (= how many weeks have passed since `min_week`).<br>\n<font color=\"RoyalBlue\">ここから、base_FVC（＝Patient の min_week 時の FVC）と base_week（＝min_week から何週経ったときのデータか）を算出していきます。</font><br>\n\n\nFirst, extract the rows in the data where `Weeks` is `min_week` and set them to `base`.\nExtract only the `Patient` and `FVC` columns from the `base`, and change the column name from `FVC` to `base_FVC`.\nThen create a new `nb` column and set all the values to 1.\nGroup the `base` with `Patient` and compute the cumulative sum with the `nb` column.\nExtract only the rows from `base` that have `nb` columns of 1, and replace `base`.\nThis allows us to eliminate duplicate `Patient` rows from the `base` dataframe with `base_FVC` in it.\nLet's remove the `nb` column.<br>\n<font color=\"RoyalBlue\">まず、data の Weeks が min_week である行を抽出し、base とします。<br>\nbase から Patient, FVC 列だけを抜き出します。<br>\n列名を FVC から base_FVC に変更します。<br>\nそして新たに nb 列を作り、値をすべて1とします。<br>\nbase を Patient でグループ化し、nb 列を指定して累積和を計算します。<br>\nbase から nb 列が1の行のみを抽出し、base を置き換えます。<br>\nこれにより、base_FVC が載っているデータフレーム base から Patient の重複を無くすことができます。<br>\nnb 列は削除しておきましょう。</font>"},{"metadata":{"trusted":true},"cell_type":"code","source":"base = data.loc[data.Weeks == data.min_week]\nbase = base[['Patient', 'FVC']].copy()\nbase.columns = ['Patient', 'base_FVC']\nbase['nb'] = 1\nbase['nb'] = base.groupby('Patient')['nb'].transform('cumsum')\nbase = base[base.nb==1]\nbase.drop('nb', axis=1, inplace=True)\n\nbase.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Next, we associate `base` with `Patient` in `data`.\nWe will create `base_week` column in `data`, which will be `Weeks` minus `min_week`.\nThis will add the `base_FVC` and `base_week` columns to `data`.\nWe should remove the `base` column.<br>\n<font color=\"RoyalBlue\">次に、data に base を Patient に紐づけて結合します。<br>\ndata に base_week 列を作成し、Weeks から min_week を引いた値とします。<br>\nこれで、data に base_FVC 列と base_week 列が追加されます。<br>\nbase は削除しておきましょう。</font>"},{"metadata":{"trusted":true},"cell_type":"code","source":"data = data.merge(base, on='Patient', how='left')\ndata['base_week'] = data['Weeks'] - data['min_week']\ndel base\n\ndata.head(10)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Perform one-hot-encoding of `Sex` and `SmokingStatus`.<br>\n<font color=\"RoyalBlue\">Sex と SmokingStatus のワンホットエンコーディングを行います。</font>"},{"metadata":{"trusted":true},"cell_type":"code","source":"categorical_features = ['Sex', 'SmokingStatus']\nfeatures_nn = []\nfor col in categorical_features:\n    for mod in data[col].unique():\n        features_nn.append(mod)\n        data[mod] = (data[col] == mod).astype(int)\n\ndata.head(10)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Normalize `Percent`, `Age`, `base_FVC`, and `base_week`.<br>\n<font color=\"RoyalBlue\">Percent, Age, base_FVC, and base_week の正規化を行います。</font>"},{"metadata":{"trusted":true},"cell_type":"code","source":"data['Percent_n'] = (data['Percent'] - data['Percent'].min() ) / ( data['Percent'].max() - data['Percent'].min() )\ndata['Age_n'] = (data['Age'] - data['Age'].min() ) / ( data['Age'].max() - data['Age'].min() )\ndata['base_FVC_n'] = (data['base_FVC'] - data['base_FVC'].min() ) / ( data['base_FVC'].max() - data['base_FVC'].min() )\ndata['base_week_n'] = (data['base_week'] - data['base_week'].min() ) / ( data['base_week'].max() - data['base_week'].min() )\nfeatures_nn += ['Age_n', 'Percent_n', 'base_week_n', 'base_FVC_n']\n\nprint(features_nn)\ndata.head(10)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now that the process is done, let's split data into `train`, `test` and `sub` using the `WHERE` column, and remove `data`.<br>\n<font color=\"RoyalBlue\">処理が終わったので、WHERE 列を使って data を train, test, sub に分割し直しましょう。data は削除しておきます。</font>"},{"metadata":{"trusted":true},"cell_type":"code","source":"train = data.loc[data.WHERE=='train']\ntest = data.loc[data.WHERE=='val']\nsub = data.loc[data.WHERE=='test']\ndel data\n\ntrain.shape, test.shape, sub.shape","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"<a id=\"quantile_m\"></a>\n## Build the model 🧠\n\nThis competition is evaluated on a modified version of the Laplace Log Likelihood. For each true FVC measurement, you will predict both an FVC and a confidence measure (standard deviation σ). The metric is computed as:<br>\n<font color=\"RoyalBlue\">このコンペでは、ラプラス対数尤度の修正版で評価されます。真の各 FVC 測定について、FVC と信頼度（標準偏差 σ）の両方を予測します。メトリックは次のように計算されます。</font><br><br>\n\n$\\large \\sigma_{clipped} = max(\\sigma, 70),$<br>\n$\\large \\Delta = min ( |FVC_{true} - FVC_{predicted}|, 1000 ),$<br>\n$\\Large metric = -   \\frac{\\sqrt{2} \\Delta}{\\sigma_{clipped}} - \\ln ( \\sqrt{2} \\sigma_{clipped} ).$<br>\n\nIn the following code, C1 is the value of confidence clipping in the modified Laplace Log Likelihood, an evaluation metric, and C2 is the error threshold.<br>\n<font color=\"RoyalBlue\">下記のコードにおいて、C1 は評価指標である修正版ラプラス対数尤度における信頼度のクリッピングの値、C2 は誤差の閾値です。</font><br>\n\nThe `score` function takes the true and predicted values of the target variable and returns a score based on the modified Laplace Log Likelihood.\nThe `qloss` function is a pinball loss function, which is the loss function used when a multiple-quantile regression prediction is trained.\nThe `mloss` function takes a percentage and returns a function that sums the return values of the `score` function and the `qloss` function according to the percentage.<br>\n<font color=\"RoyalBlue\">score 関数はこの評価メトリックです。コンペの目的変数の真の値と予測値を受け取り、修正版ラプラス対数尤度に基づいたスコアを返します。<br>\nqloss 関数は、重分位点回帰予測が学習するときに使用する損失関数であるピンボールロス関数です。<br>\nmloss 関数は割合 _lambda を受け取り、その割合に応じて score 関数と qloss 関数の戻り値を合計する関数を返します。</font><br><br>\nHere, we define confidence as the difference between the predicted values at 0.2 and 0.8 quartiles.<br>\n<font color=\"RoyalBlue\">ここで、信頼度を0.2分位点と0.8分位点における予測値の差として定義しています。</font>"},{"metadata":{"trusted":true},"cell_type":"code","source":"C1, C2 = tf.constant(70, dtype=\"float32\"), tf.constant(1000, dtype=\"float32\")\n#=============================#\ndef score(y_true, y_pred):\n    tf.dtypes.cast(y_true, tf.float32)\n    tf.dtypes.cast(y_pred, tf.float32)\n    sigma = y_pred[:, 2] - y_pred[:, 0]\n    fvc_pred = y_pred[:, 1]\n\n    sigma_clip = tf.maximum(sigma, C1)\n    delta = tf.abs(y_true[:, 0] - fvc_pred)\n    delta = tf.minimum(delta, C2)\n    sq2 = tf.sqrt( tf.dtypes.cast(2, dtype=tf.float32) )\n    metric = (delta / sigma_clip)*sq2 + tf.math.log(sigma_clip * sq2)\n    return backend.mean(metric)\n#============================#\ndef qloss(y_true, y_pred):\n    # Pinball loss for multiple quantiles\n    qs = [0.2, 0.5, 0.8]\n    q = tf.constant(np.array([qs]), dtype=tf.float32)\n    e = y_true - y_pred\n    v = tf.maximum(q*e, (q-1)*e)\n    return backend.mean(v)\n#=============================#\ndef mloss(_lambda):\n    def loss(y_true, y_pred):\n        return _lambda * qloss(y_true, y_pred) + (1 - _lambda)*score(y_true, y_pred)\n    return loss","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def make_model():\n    inp = L.Input(len(features_nn), name=\"Patient\")\n    x = L.Dense(100, activation=\"relu\", name=\"d1\")(inp)\n    x = L.Dense(100, activation=\"relu\", name=\"d2\")(x)\n    p1 = L.Dense(3, activation=\"linear\", name=\"p1\")(x)\n    p2 = L.Dense(3, activation=\"relu\", name=\"p2\")(x)\n    preds = L.Lambda(lambda x: x[0] + tf.cumsum(x[1], axis=1), \n                     name=\"preds\")([p1, p2])\n    \n    model = M.Model(inp, preds, name=\"NeuralNet\")\n    model.compile(loss=mloss(0.64),    # changed from 0.8\n                  optimizer=tf.keras.optimizers.Adam(lr=0.1, decay=0.01),\n                  metrics=[score])\n    return model\n\nmodel = make_model()\nmodel.summary()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_model(model)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"<a id=\"quantile_c\"></a>\n## Cross validation 💭\n\nUse the model we created to cross-validate.<br>\n<font color=\"RoyalBlue\">作成したモデルを使ってクロスバリデーションを行いましょう。</font>"},{"metadata":{"trusted":true},"cell_type":"code","source":"X_train = train[features_nn].values\nX_test = sub[features_nn].values\n\ny_train = train['FVC'].values\n\noof_train = np.zeros((X_train.shape[0], 3))\ny_preds = np.zeros((X_test.shape[0], 3))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"BATCH_SIZE = 128\nEPOCHS = 804    # changed from 800\nNFOLD = 5\n\nkf = KFold(n_splits=NFOLD)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\nfor fold_id, (tr_idx, va_idx) in enumerate(kf.split(X_train)):\n    print(f\"FOLD {fold_id+1}\")\n    model = make_model()\n    model.fit(X_train[tr_idx], y_train[tr_idx], batch_size=BATCH_SIZE, epochs=EPOCHS, \n              validation_data=(X_train[va_idx], y_train[va_idx]), verbose=0)\n    print(\"train\", model.evaluate(X_train[tr_idx], y_train[tr_idx], verbose=0, batch_size=BATCH_SIZE))\n    print(\"val\", model.evaluate(X_train[va_idx], y_train[va_idx], verbose=0, batch_size=BATCH_SIZE))\n    oof_train[va_idx] = model.predict(X_train[va_idx], batch_size=BATCH_SIZE, verbose=0)\n    y_preds += model.predict(X_test, batch_size=BATCH_SIZE, verbose=0) / NFOLD","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's illustrate the correct and predicted values.<br>\n<font color=\"RoyalBlue\">正解値と予測値を図示してみましょう。</font>"},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(figsize=(12, 12))\n\nidxs = np.random.randint(0, y_train.shape[0], 100)\nax.plot(y_train[idxs], label=\"ground truth\", color=palette_ro[0])\nax.plot(oof_train[idxs, 0], label=\"q20\", color=palette_ro[3], ls=':', alpha=0.5)\nax.plot(oof_train[idxs, 1], label=\"q50\", color=palette_ro[4], ls=':', alpha=0.5)\nax.plot(oof_train[idxs, 2], label=\"q80\", color=palette_ro[5], ls=':', alpha=0.5)\nax.legend(loc=\"best\");","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We calculate the optimized 𝜎 (standard deviation) from the `oof_train`. `sigma_opt` is the mean absolute error between the correct value of each fold and the prediction (median), `sigma_unc` is the difference between the prediction (0.2 quantile) and the prediction (0.8 quantile), and `sigma_mean` is the mean value of the difference.<br>\n<font color=\"RoyalBlue\">では、oof_train から最適化された 𝜎（標準偏差）を計算しましょう。各フォールドの正解値と予測値（中央値）との平均絶対誤差を sigma_opt, 予測値（0.2分位数）と予測値（0.8分位数）との差を sigma_unc, その平均値を sigma_mean とします。</font>"},{"metadata":{"trusted":true},"cell_type":"code","source":"sigma_opt = mean_absolute_error(y_train, oof_train[:, 1])\nsigma_unc = oof_train[:, 2] - oof_train[:, 0]\nsigma_mean = np.mean(sigma_unc)\nprint(sigma_opt, sigma_mean)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(sigma_unc.min(), sigma_unc.mean(), sigma_unc.max(), (sigma_unc>=0).mean())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(np.mean(y_train / oof_train[:, 1]))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(figsize=(16, 6))\n\nsns.distplot(sigma_unc, ax=ax, color=palette_ro[1])\nax.set_title(\"uncertainty in prediction\", fontsize=18);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"<a id=\"submit\"></a>\n# Ensemble & Submit 📝"},{"metadata":{"trusted":true},"cell_type":"code","source":"sub.head(10)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Prepare a submission file from the predictions of the neural network.<br>\n<font color=\"RoyalBlue\">ニューラルネットワークの予測結果から提出ファイルを作成する準備をします。</font>"},{"metadata":{"trusted":true},"cell_type":"code","source":"sub['FVC1'] = y_preds[:, 1]\nsub['Confidence1'] = y_preds[:, 2] - y_preds[:, 0]\n\nsub.head(10)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"`subm` is defined by extracting the required columns from the `sub` and leaving only the rows with non-null data in `FVC1`.<br>\n<font color=\"RoyalBlue\">`sub` から必要な列を抜き出し、`FVC1` のデータが `null` でない行だけにしたものを `subm` とします。</font>"},{"metadata":{"trusted":true},"cell_type":"code","source":"subm = sub[['Patient_Week', 'FVC', 'Confidence', 'FVC1', 'Confidence1']].copy()\nsubm.loc[~subm.FVC1.isnull(),'FVC'] = subm.loc[~subm.FVC1.isnull(),'FVC1']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"if sigma_mean<70:\n    subm['Confidence'] = sigma_opt\nelse:\n    subm.loc[~subm.FVC1.isnull(),'Confidence'] = subm.loc[~subm.FVC1.isnull(),'Confidence1']\n\nsubm.head(10)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"subm.describe().T","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Read the original `test.csv` and overwrite the `FVC` and `Confidence` in the predicted data to be submitted if they are known in the `test.csv`.<br>\n<font color=\"RoyalBlue\">オリジナルの test.csv を読み込み、投稿予定の予測データの中に test.csv で既知のデータがあれば FVC と Confidence を上書きします。</font>"},{"metadata":{"trusted":true},"cell_type":"code","source":"org_test = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/test.csv')\nfor i in range(len(org_test)):\n    subm.loc[subm['Patient_Week']==org_test.Patient[i]+'_'+str(org_test.Weeks[i]), 'FVC'] = org_test.FVC[i]\n    subm.loc[subm['Patient_Week']==org_test.Patient[i]+'_'+str(org_test.Weeks[i]), 'Confidence'] = 70\n\nsubm[[\"Patient_Week\",\"FVC\",\"Confidence\"]].to_csv(\"submission_regression.csv\", index=False)\nreg_sub = subm[[\"Patient_Week\",\"FVC\",\"Confidence\"]].copy()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Ensemble two models.<br>\n<font color=\"RoyalBlue\">２つのモデルをアンサンブルします。</font>"},{"metadata":{"trusted":true},"cell_type":"code","source":"df1 = img_sub.sort_values(by=['Patient_Week'], ascending=True).reset_index(drop=True)\ndf2 = reg_sub.sort_values(by=['Patient_Week'], ascending=True).reset_index(drop=True)\n\ndf = df1[['Patient_Week']].copy()\ndf['FVC'] = 0.2*df1['FVC'] + 0.8*df2['FVC']    # changed from 0.25, 0.75\ndf['Confidence'] = 0.0*df1['Confidence'] + 1.0*df2['Confidence']    # changed from 0.26, 0.74\ndf.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df.to_csv('submission.csv', index=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We could create `submission.csv`. Thank you so much for reading!<br>\n<font color=\"RoyalBlue\">submission.csv を作成できました。読んでくださりありがとうございました！</font>"}],"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}