{"cells":[{"metadata":{},"cell_type":"markdown","source":"## Basic imports"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport os\nimport random\nfrom tqdm.auto import tqdm\n\n\nimport matplotlib.pyplot as plt\nimport cv2\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.metrics import mean_absolute_error\nfrom sklearn.model_selection import KFold, RepeatedKFold\n\nimport tensorflow as tf\nimport tensorflow.keras.backend as K\nimport tensorflow.keras.layers as L\nimport tensorflow.keras.models as M","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"seed = 123\n\ndef seed_everything(seed=123):\n random.seed(seed)\n os.environ['PYTHONHASHSEED'] = str(seed)\n np.random.seed(seed)\n tf.random.set_seed(seed)\n \nseed_everything(seed)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# data params\nmean = -0.13859619186395233 #from training notebook\nstd = 0.1546567782480763 #from training notebook\nnum_cv_folds = 4\nnum_repeats = 50","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Quantile Regression\nSee: https://www.kaggle.com/ulrich07/osic-multiple-quantile-regression-starter"},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Data preparation**\n\nSee: https://www.kaggle.com/mattbast/feature-engineering-with-a-linear-model"},{"metadata":{"trusted":true},"cell_type":"code","source":"## construct train input\ntrain_data = pd.read_csv(\"../input/osic-pulmonary-fibrosis-progression/train.csv\")\nprint(train_data.shape)\n\n# merge \"base values\"\nbase_df = train_data.copy()[['Patient', 'FVC', 'Percent', 'Weeks']].sort_values(['Patient','Weeks']).groupby('Patient').head(1)\nbase_df.columns = ['Patient', 'base_FVC', 'base_Percent', 'base_Weeks']\ntrain_data = train_data.merge(base_df, on='Patient')\nprint(train_data.shape)\n\n\ntrain_data['Weeks_passed'] = train_data.Weeks-train_data.base_Weeks\ntrain_data = train_data[train_data.Weeks_passed>0]\n#train_data = train_data.groupby(['Patient','base_Weeks']).tail(10)\ntrain_data = train_data.drop_duplicates()\nprint(train_data.shape)\n\n\n# folds = pd.read_csv('../input/osic-targets/train_folds.csv')[['Patient', 'split-all']].drop_duplicates()\n# folds.columns = ['Patient', 'fold']\n# folds['fold'] = folds['fold'].astype(int)\n# train_data = train_data.merge(folds)\n# print(train_data.shape)\n\ntrain_data['est_Percent'] = train_data.apply(lambda x: mean*(x.Weeks-x.base_Weeks) + x.base_Percent, axis=1)\ntrain_data['est_FVC'] = train_data['base_FVC']/train_data['base_Percent']*train_data['est_Percent']\n\ntrain_data = train_data[['Patient', 'Weeks', 'FVC', 'Percent', 'Age', 'Sex', 'SmokingStatus', \n 'base_FVC', 'base_Percent', 'base_Weeks', 'Weeks_passed', #'fold',\n 'est_Percent', 'est_FVC',\n ]].drop_duplicates()\ntrain_data['WHERE'] = 'train'\n\nprint(train_data.shape)\ntrain_data.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_data.describe().T","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sample_sub = pd.read_csv(\"../input/osic-pulmonary-fibrosis-progression/sample_submission.csv\")\nsample_sub['Patient'] = sample_sub['Patient_Week'].apply(lambda x: x.split(\"_\")[0])\nsample_sub['Weeks'] = sample_sub['Patient_Week'].apply(lambda x: x.split(\"_\")[1])\nsample_sub.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# construct test_data input\ntest_data = pd.read_csv('../input/osic-pulmonary-fibrosis-progression/test.csv')\ntest_data = test_data.groupby('Patient').first()[['FVC', 'Percent', 'Weeks', 'Percent', 'Age', 'Sex', 'SmokingStatus',]]\ntest_data.columns = ['base_FVC', 'base_Percent', 'base_Weeks', 'Percent', 'Age', 'Sex', 'SmokingStatus',]\n\ntest_data = test_data.merge(sample_sub[['Patient','Weeks']], on='Patient')\ntest_data['Weeks'] = test_data.Weeks.astype(int)\ntest_data['Weeks_passed'] = test_data.Weeks-test_data.base_Weeks\n\nprint(test_data.shape)\ntest_data['est_Percent'] = test_data.apply(lambda x: mean*(x.Weeks-x.base_Weeks) + x.base_Percent, axis=1)\ntest_data['est_FVC'] = test_data['base_FVC']/test_data['base_Percent']*test_data['est_Percent']\n \nprint(test_data.shape)\n \ntest_data = test_data[['Patient', 'Weeks', 'Age', 'Sex', 'SmokingStatus', \n 'base_FVC', 'base_Percent', 'base_Weeks', 'Weeks_passed', \n 'est_Percent', 'est_FVC',]].drop_duplicates()\ntest_data['WHERE'] = 'test'\nprint(test_data.shape)\ntest_data.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_data.describe().T","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Combine datasets**"},{"metadata":{"trusted":true},"cell_type":"code","source":"def calculate_height(row):\n if row['Sex'] == 'Male':\n return row['base_FVC'] / (27.63 - 0.112 * row['Age'])\n else:\n return row['base_FVC'] / (21.78 - 0.101 * row['Age'])\n \ndef scale_feature(series):\n return (series - series.min()) / (series.max() - series.min())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"all_data = train_data.append(test_data)\nall_data['Height'] = all_data.apply(calculate_height, axis=1)\n\n# to categorical\nall_data = pd.concat([\n all_data,\n pd.get_dummies(all_data.Sex),\n pd.get_dummies(all_data.SmokingStatus)\n], axis=1)\n\nall_data = all_data.drop(columns=['Sex', 'SmokingStatus'])\n \n\nall_data.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# normalize\nall_data['normalized_Weeks'] = scale_feature(all_data['Weeks'])\nall_data['normalized_base_Percent'] = scale_feature(all_data['base_Percent'])\n\nall_data['normalized_Age'] = scale_feature(all_data['Age'])\nall_data['normalized_base_Weeks'] = scale_feature(all_data['base_Weeks'])\nall_data['normalized_base_FVC'] = scale_feature(all_data['base_FVC'])\nall_data['normalized_Weeks_passed'] = scale_feature(all_data['Weeks_passed'])\nall_data['normalized_Height'] = scale_feature(all_data['Height'])\nall_data['normalized_est_FVC'] = scale_feature(all_data['est_FVC'])\nall_data['normalized_est_Percent'] = scale_feature(all_data['est_Percent'])\n\nall_data.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"FEATURE_COLS = ['Female', 'Male', \n 'Currently smokes', 'Ex-smoker', 'Never smoked',\n 'normalized_Weeks', 'normalized_base_Weeks', 'normalized_Weeks_passed', \n 'normalized_base_Percent',\n 'normalized_base_FVC', \n 'normalized_Age', \n #'normalized_Height',\n # 'normalized_est_FVC', 'normalized_est_Percent',\n ]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"tr = all_data.loc[all_data.WHERE=='train'].reset_index(drop=True)\nte = all_data.loc[all_data.WHERE=='test'].reset_index(drop=True)\n\ntr.shape, te.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"C1, C2 = tf.constant(70, dtype='float32'), tf.constant(1000, dtype=\"float32\")\n#=============================#\ndef score(y_true, y_pred):\n tf.dtypes.cast(y_true, tf.float32)\n tf.dtypes.cast(y_pred, tf.float32)\n sigma = y_pred[:, 2] - y_pred[:, 0]\n fvc_pred = y_pred[:, 1]\n \n #sigma_clip = sigma + C1\n sigma_clip = tf.maximum(sigma, C1)\n delta = tf.abs(y_true[:, 0] - fvc_pred)\n delta = tf.minimum(delta, C2)\n sq2 = tf.sqrt( tf.dtypes.cast(2, dtype=tf.float32) )\n metric = (delta / sigma_clip)*sq2 + tf.math.log(sigma_clip* sq2)\n return -K.mean(metric)\n\ndef qloss(y_true, y_pred):\n # Pinball loss for multiple quantiles\n qs = [0.2, 0.50, 0.8]\n q = tf.constant(np.array([qs]), dtype=tf.float32)\n e = y_true - y_pred\n v = tf.maximum(q*e, (q-1)*e)\n return K.mean(v)\n\ndef mloss(_lambda):\n def loss(y_true, y_pred):\n return _lambda * qloss(y_true, y_pred) + (1 - _lambda)*score(y_true, y_pred)\n return loss\n\n#=================\ndef make_model(nh):\n z = L.Input((nh,), name=\"Patient\")\n x = L.Dense(100, activation=\"relu\", name=\"d1\")(z)\n x = L.Dense(100, activation=\"relu\", name=\"d2\")(x)\n p1 = L.Dense(3, activation=\"linear\", name=\"p1\")(x)\n p2 = L.Dense(3, activation=\"relu\", name=\"p2\")(x)\n preds = L.Lambda(lambda x: x[0] + tf.cumsum(x[1], axis=1), \n name=\"preds\")([p1, p2])\n \n model = M.Model(z, preds, name=\"CNN\")\n model.compile(loss=mloss(0.65), optimizer=tf.keras.optimizers.Adam(lr=0.1, beta_1=0.9, beta_2=0.999, epsilon=None, decay=0.01, amsgrad=False), metrics=[score])\n return model","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"nh = len(FEATURE_COLS)\nnet = make_model(nh)\nprint(net.summary())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"NOISE_ON_FEATURES = ['FVC',\n 'normalized_Weeks', 'normalized_base_Weeks', 'normalized_Weeks_passed',\n 'normalized_base_Percent', 'normalized_base_FVC', \n 'normalized_Age', 'normalized_Height',\n #'normalized_est_FVC', 'normalized_est_Percent'\n ]\n\nNOISE_MAX_PCT = 0.05\nnoise_generator = lambda x: x*(1+(np.random.randint(-NOISE_MAX_PCT*1000, NOISE_MAX_PCT*1000)/1000))\n\nN_ARTIFICIAL_SAMPLES = len(tr)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"BATCH_SIZE = 128\nEPOCHS = 200","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"pred = np.zeros((len(tr), 3))\npe = np.zeros((len(te[FEATURE_COLS]), 3))\n\n\nrkf = RepeatedKFold(n_splits=num_cv_folds, n_repeats=num_repeats, random_state=seed)\nfold = 0\nfor tr_idx, val_idx in rkf.split(tr):\n print(f\"FOLD {fold}\")\n X_tr = tr.loc[tr_idx]\n X_tr_artif = X_tr.sample(n=N_ARTIFICIAL_SAMPLES, replace=True)\n X_tr_artif[NOISE_ON_FEATURES] = X_tr_artif[NOISE_ON_FEATURES].apply(noise_generator)\n X_tr = X_tr.append(X_tr_artif)\n y_tr = X_tr['FVC'].values\n X_tr = X_tr[FEATURE_COLS].values\n \n X_val = tr.loc[val_idx]\n #X_val[NOISE_ON_FEATURES] = X_val[NOISE_ON_FEATURES].apply(noise_generator)\n y_val = X_val['FVC'].values\n X_val = X_val[FEATURE_COLS].values\n \n \n X_te = te[FEATURE_COLS].values\n \n net = make_model(nh)\n net.fit(X_tr, y_tr, \n batch_size=BATCH_SIZE, epochs=EPOCHS, \n validation_data=(X_val, y_val), verbose=0) #\n \n print(\"train\", net.evaluate(X_tr, y_tr, verbose=0, batch_size=BATCH_SIZE))\n print(\"val\", net.evaluate(X_val, y_val, verbose=0, batch_size=BATCH_SIZE))\n print(\"predict val...\", end=\" \")\n pred[val_idx] += net.predict(X_val, batch_size=BATCH_SIZE, verbose=0)\n print(\"done\")\n \n print(\"predict test...\", end=\" \")\n pe += net.predict(X_te, batch_size=BATCH_SIZE, verbose=0) / (num_cv_folds*num_repeats)\n print(\"done\")\n fold+=1\n#==============\n\npred = pred/num_repeats","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X = tr[FEATURE_COLS].values\ny = tr['FVC'].values\n\nsigma_opt = mean_absolute_error(y, pred[:, 1])\nunc = pred[:,2] - pred[:, 0]\nsigma_mean = np.mean(unc)\nprint(sigma_opt, sigma_mean)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"tr['FVC_pred'] = pred[:, 1]\ntr['Confidence'] = pred[:,2] - pred[:, 0]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def lb_metric(train):\n train['sigma_clipped'] = train['Confidence'].apply(lambda x: max(x, 70))\n train['diff'] = abs(train['FVC'] - train['FVC_pred'])\n train['delta'] = train['diff'].apply(lambda x: min(x, 1000))\n train['score'] = -np.sqrt(2)*train['delta']/train['sigma_clipped'] - np.log(np.sqrt(2)*train['sigma_clipped'])\n score = train['score'].mean()\n return score\n\n\nscore = lb_metric(tr)\nprint(f'Local Score: {score}')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"idxs = np.random.randint(0, y.shape[0], 100)\nplt.figure(figsize=(16,9))\nplt.plot(y[idxs], label=\"ground truth\")\nplt.plot(pred[idxs, 0], label=\"q25\")\nplt.plot(pred[idxs, 1], label=\"q50\")\nplt.plot(pred[idxs, 2], label=\"q75\")\nplt.legend(loc=\"best\")\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(unc.min(), unc.mean(), unc.max(), (unc>=0).mean())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.hist(unc)\nplt.title(\"Uncertainty in prediction\")\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(5, 1, figsize=(10, 20))\n\nfor i, pt in enumerate(np.random.choice(tr['Patient'], 5)):\n patient_log = tr[tr['Patient'] == pt]\n\n ax[i].set_title(pt)\n ax[i].scatter(patient_log['Weeks_passed'], patient_log['FVC'], label='truth')\n ax[i].scatter(patient_log['Weeks_passed'], patient_log['FVC_pred'], label='prediction')\n ax[i].legend()\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"xx = np.linspace(70, 1000, 100)\nbest_y = -12\n\nfor x in xx:\n tr['Confidence'] = x\n yy = lb_metric(tr)\n plt.scatter(x, yy, color='k')\n \n if best_y